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ABSTRACT 


This investigation explores the effects of chordwise forces and 
deformations and steady-state deformation due to lift on the static 
and dynamic aeroelastic stability of a uniform cantilever wing. Results 
of this analysis are believed to have practical applications for high- 
performance sailplanes and certain RPV's- 

The airfoil cross section is assumed to be symmetric nnd camber 
bending is neglected. Motions in vertical bending, fore-and-aft 
bending, ana torsion are considered. A differential equation model is 
developed, which includes the nonlinear elastic bending- tors ion coupling 
that accompanies even moderate deflections. A linearized expansion in 
small time-dependent perturbation deflections is made about a steady 
flight condition. The stability determinant of the linearized system 
then contains coefficients that depend on steady displacements. Loads 
derived from two-dimensional incompressible aerodynamic theory are used 
to obtain the majority of the results, but cases using three-dimensional 
subsonic compressible theory are also studied. 

The stability analysis is carried out in terms of the dynamically 
uncoupled natural modes of vibration of the uniform cantilever. Dynamic 
stability in the case of incompressible strip-theory airloads is deter- 
mined in two ways. One is the "V-g method" familiar to aeroelastlcians. 


When steady deformations are present this method requires an iterative ' " '' 
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matching of flutter speeds with estimated speed. The second approach ' O 


involves determination of the complex eigenvalues of the aeroelastic 
modes at any desired flight condition. The aerodynamic loads are 
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expressed Ju terms of the {generalized Theodorsen function; eif*envalues 
of the aeroelastii system are located with a gradient search technique. 

The effect of steady drag on flutter of nonlifting wings using 
IncomFireasibie strip-theory Is studied and shown to correlate with 
previously known results. Next, the Influence of steady lifting 
deformations on flutter Is Investigated, and flutter modes are found 
that involve fore-and-aft bending motions. The significance of unsteady 
leading edge suction forces, which are predicted by the two-dimensional 
incompressible aerodynamic theory. Is then examined. Two idealized 
examples based upon existing sailplanes are analyzed. 

Steady drag loads lower the flutter speed for larger aspect ratios 
hut increase it for aspect rattca less than a certain value. Divergence 
speed is more sensitive to steady drag, and for very high aspect ratio 
wings it can fall below the bending-torsion flutter speed. Steady 
deformations due to lift always decrease the flutter speed by an amount 
dependent upon the aspect ratio and the fore-and-aft bending stiffness. 
Leading-edge suction forces increase Clutter speed. 

Three-dimensional steady and un.steady airloads are introduced into 
the V-g flutter analysis scheme, and for a ^fach number of zero the role 
of steady lifting deformations and unsteady leading-edge suction forces 
Is more accurately determined.. The behavior predicted using strip 
theory loads is again observed, and the suction forces are confirmed to 
contribute a significant stabilizing effect. Further calculations using 
high subsonic Mach numbers reveal only mild effects due to compressi- 
bility (disregarding unsteady chordwise loads of viscous origin). 
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Chapter 1 


INTRODUCTION 

Aerodynamic loads and deformations parallel to the chord ate 
usually neRlected during studies of aeroelastic stability of lifting 
surfaces. Furthermore, dynamic stability is usually analyzed 'Without 
regard for the influence of steady deformations due to steady -state 
lift. 

Very little literature exists »/hich treats chordwise f'/»rces and 
bending in an a-rtroelastif.- analysis. A substantial study oH the effects 
of drag lunus OR diverg'ince of a cantilever wing Is made ay Petre (R/af. 
2, pp. 449-487). Here it is clearly der^ionscraced that ' he Interaction 
of drag with bending deformations due to lifting loads can signifi- 
cantly reduce diverr/ence speeds. Goetz (Ref. 17) co/isidered this same 
drag-bending deformation divergence mechanism, specialized to the case 
of a rigid lifting surface at the end of a beam-rod. This work involved 
supersonic flow past a surface having a blunt leading edge, and the 
resulting sizeable drag forces caused a significant reduction of the 
classical divergence speed. 

One example of an aeroelastic study in which chordwise deformations 
of a straight cantilever wing are accounted for is the work in the area 
of tilting proprotor aircraft by Wayne Johnson. The cruising flight 
condition (Ref. 27) is modeled using a proprotor with axial flow mounted 
at the tip of a cantilever wing. The additional degrees of freedom 
associated with the individual elastic rotor blades and the aerodynamic 
and inertial effects: of the proprotor result in a much more complex and 
specialized aeroelastic analysis chan Is considered here. 
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The effect which steady chordwise forces can have upon dynamic 
stability was explored by Petre (Ref. 3) and Petre and Ashley (Ref. 1) 
using two-dimensional incompressible unsteady lifting airloads. The 
latter work presents extensive calculation^^ regarding the effect of 
steady drag on bending-torsion flutter of a uniform cantilever wing. 

It serves aa the starting point for the work pursued in this thesis. 

The objectives of this thesis are as follows: 

1. To check .nnd interpret the predicted effect of steady drag 
on the flutter behavior of a nonllfting wing discovered in 
Ref. 1, using a modal approach Instead of a collocation 
.approach. 

2. To generalize the equations of motion to include consistently 
fore-and-aft bending motions, adequately accounting for the 
elastic coupling among the three degrees of freedom. 

3. To Include steady-state lifting deformations in the dynamic 
stability analysis by considering small time-dependent per- 
turbation deflections about a steady displacement solution. 

4. To allow for unsteady leading-edge-suction forces in the 
chordwise direction predicted by two-dimensional incompressible 
unsteady potential flow theory. 

5. To Improve the representation of both steady and unsteady air- 
loads by use of a three-dimensional subsonic kernel function 
program, from which leading-edge suction and induced drag can 
also be obtained. 

Items 2, 3, and 4 arc interrelated and together represent a consis- 
tent extension over previous research in the modeling of the physical 
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svsccB. For example.', when sceadv lateral bending detoraat U>n!i are 
present, unsceadv chordwise loads can induce cwiKtlng stot ions which 
slgnif leant Iv affect dvnaniic stabilltv. The naln purpose here will 
be CO establish trends and gain fundanental Insight into the influence 
of cherdwise forces and sCcadv detonsutt Ions, hopeful Iv shedding light 
on their importance in priictical aerospace problems. 

Certain assumptions are adhered to throughout this thesis. Since 
the emphasis is upon working from the equ.itions of siot ion in differen- 
tial form in order to include certain nonlinear elastic coupling terms, 
it is convenient to restrict this studv to straight cantilever wings 
having mass and stiffness properties uniform with span. Tlio wing is 
taken to be a one-dimensional structure in the sense that all dofonaa- 
tions are described as functions of the spanwtse v.ariahle v . Camber 
bending is neglected and the simple Euler-Bemoul 1 I beam stress-strain 
assumptions itre used. The platform is rectangular, and the steady and 
unsteady flow fields are always assumed to be superposnble ; vmsteadv 
loads are computed for the undefonned geometry and applied to the deformed 
wing. Although these assumptions would be restrictive for the purpose 
of modeling actual structures, they are acceptable here since only the 
relative influences of chordwise loads and steady deformations are of 
interest . 

In Chapter II a vertical-bending/torsion model basic to the system 
of Ref. 1 is developed and modal equations are derived to permit flutter 
calculations for zero steady .lift with steady drag included. Assumed 
mode solutions are then compared with results of Ref. 1. A linear 
steady-state version of the modal equations is then examined to allow 
computation of divergence speeds as affected by steady drag. 
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In Chapter III a model central to this thesis Is developed that 
Includes vertical bending, chordwise bending, and torsion. Nonlinear 
steady and linearized unsteady differential equations are then deduced, 
and they are analytically compared with the model in Chapter II. 

In Chapter IV the modal forms of these steady and linearized 
unsteady equations are set up to include lifting airloads deduced from 
Incompressible steady and unsteady strip-theory, k scheme based on the 
scv-called V-g nwithod of flutter analysis is used to determine neutral 
dynamic stability conditions, and results are checked against those of 
Chapter II. 

In Chapter V a generalization to the case of arbitrary motion in 
tine Is presented through Laplace transformation of the modal equations, 
which requires Incompressible unsteady two-dimensional airloads valid 
for non-perlodlc motions of the wing. 

A determinant iteration procedure Is used to determine the aero- 
elastic eigenvalues for flight speeds .ibove and below the flutter speed. 
Fir, ally, the effect of unsteady leading-edge auction forces predicted 
by inconpresolble strip theory is included In the linearized unsteady 
stability system. 

In Chapter VI all results for incompressible strip-theory airloads 
are assembled and systematically presented, concluding with tb#o examples 
based upon actual sailplanes. 

In Chapter VII the flutter speed prediction scheme of Chapter TV 
is modified to use three-dimensional subsonic steady and unsteady ' 

ll,fting airloads. Results are presented to indicate the effects of 
three-dimensional aerodynamics, unsteady drag, and compressibility. 
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In Chapter VIII the principal conclusions of this thesis are 
presented. 
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Chapter II 


DRAG EFFKCTS ON FLUTTER AND DIVERGENCE; 
DEFORMATION IN VERTICAL BENDING AND TORSION ONLY 


A . The Vertical Bending/Torston Equations 

In the Petre and AshJ.ey work (Ref. 1) the effect of drag on Elutl-or 
of cantilever wings is studied assuming only vertical bending and 
torsion about tlie elastic axis. The equations are cast into integral 
form, and solutions are obtained by collocation of the integral equations 
at ten stations across the span. Here wc check the results of Ref. 1 
using an entirely different, modal approach. The results also servo as 
examples with which to compare solutions found with the more gener.al 
system, developed in the following chapter, which includes chordwise 
bending. 

The differential equations of bonding and torsion, as given in 
Ref. 1 and adapted to the present noLaCic.ri and coordinate system (as 
in Fig. 2-1), are 


(2-la) 


[hl^ V"’ ' ^ ^ 


+ ^ inj) = 0 


(2-lb) ^ [GI^ |i] _ j 0 + 3^ 1^ + 


M ri 

M -K— T = 0 

z 9v‘ 


^®scsiai(o bi ailjc 
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The quantity is a function of y given by 

(2-2) M = - n(n)(n-y)dn 

z y 

= - D Vs[S.-y]^ 

where it Is assumed that the running drag force D has constant magni- 
tude across the span. can be recognized as the total moment about 

the vertical axis applied at station y by all drag acting outboard 
of this station. Kean values of the twist and bending displacement 
are assumed to be zero. 

Equations (2-1) are linear and are coupled inertially by the third 
term In each; these terms arise from the offset of the center of mass 
of the airfoil section from the elastic axis. A second coupling effect 
Is due entirely to drag. Introduced through the terms containing . 

All remaining terms in the equations represent the conventional elastic, 
inertial, and unsteady lifting aerodynamic load contributions. 

The manner in which the drag loads have been introduced into the 
system is discussed by Petre (pp. A49-A87, Ref. 2) and can be explained 
by the following physically oriented argument. The drag coupling term 
in the bending equation arises from the resultant bending moment 
applied at station y due to drag forces outboard of y . As shown in 
Fig. 2-1, a drag force acting outboard of station y has a component 
D<t)(y) perpendicular to the airfoil's principal axis of vertical bending, 
giving a resultant moment at y about the principal direction of 

m (y) (J)(y)D (n-y)dn 

X y 

= (l)(y)D 

-B -M (J* 
z 


8 



The moment-curvature relation for the beam. 


-2 m 

3 w _ X 

“ El 

X 


together with moment-shear equilibrium 



= q(y) 


(here q(y) is the positive-upward running load in the principal 
direction of the section) allows the drag effect to be expressed in 
equilibrium with the elastic term as 



This is the same drag coupling term appearing in (2-1). 

A similar derivation reveals the origin c>f the drag term In the 
torsion equation. Looking at the front view of Che wing (Fig. 2.2) one 
can see that a drag load at ri gives rise to a twisting moment at y 
by acting through the moment arm given by distance e , 


e = w(ri) 



w(y) + 


3w 

3y 


(y)(h-y) 


1 

J 


The resultant torque applied at station y by all drag forces acting 
outboard is given by 

T <y) = D [w(n) - w(y) - 1^ (n-y)ldn 

y y dy 

Differentiating with respect to y , 
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3tf aw~ | 

3y 3yJ 


dn 


= - B ^ 


= M 


3*w 


z 3^ 


Elastic equilibrium for a rod loaded with applied torsion Ty(y) is 
given by 



= 0 
3y 


The drag coupling term of equation (2-1) can then be identified in 
the result 



B, Solution By Assumed Modes 

When structural dynamics problems yield solutions whose frequencies 
are within the range of the structure’s lowest normal mode natural 
frequencies and the latter have a sufficiently sparse distribution, 
modal analysis methods prove to be effective. Primary bending-torsion 
flutter of cantilever wings, under study here, is a classic example in 
aeroelasticity of such a system. Inclusion of steady drag effects 
should have only a small effect on the range of frequencies over which 
flutter solutions occur and, by this reasoning, should not adversely 
affect the convergence of modal solutions. Petre in Ref. 1 expresses 
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I 


f 


Che opinion that methods relying on the assumption of a few prescribed 
deformation modes constrain the flutter mode shape and should be avoided. 
This opinion is tested herein by actually studying modal convergence. 

Solutions of (2-1) are sought using superposition of the dynami- 
cally uncoupled natural modes of free vibration in bending and torsion 
of a uniform cantilever beam. These assumed modes, with their relevant 
properties, are described in Appendix A. Although not true normal modes 
of the inertially-coupled structure, they can he considered as "pseudo- 
orthogonal" since integrals of the type 

f f dv = 0 ± ^ i 

0 • ■' 

fh f . <r dy = 0 1 .i 

0 (Jij (Jij 

lead to uncoupled elastic behavior and hence a diagonal stiffness matrix 
in the matrix eigenvalue problem. Use of the actual normal modes would 
require that they he calculated for each wing configuration studied. 

Since the assumed modes satisfy the naLural boundary conditions at the 
free end of the cantilever as well as the geometric boundary conditions 
at its clamped root, Galerkin's method can be applied to the differential 
equations (Ref. 4, p. 218) to obtain the system in terms of modal 
generalized coordinates. 

To find the velocity for neutral stability (the flutter velocity), 
the V-g method (p. 381, Ref. 5), common in aeroelastic stability 
analyses, In employed. With simple harmonic motion of frequency (u 
the unsteady, incompressible, strip-theory lifting airloads are 
expressible as [(4-123) and (4-124) of Ref. 5] 
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(2-3a) L^(w,il);t) ■» - II p ^ + (L^ -C^a)L^14)e^'^^ 


(2-3b) tn (w,(j);t) •» H p (ij^{-(M -(‘i+a)L ] ~ 
a w vr D 


where 


+ - (V-a)(L^ + + C^a)^ 


(2-4) 


L « 1 - -^ C(k) 
w k 


S ° “ k ^ 


M - >i 
w 


M = 1 i 
^(tJ 8 k 


Here C(k) is the familiar Theodoraen function of reduced frequency 

, u)b 
“ V ■ 

The drag coupling terms in (2-1) can be treated as applied loads 
In developing the modal equations, by defining total applied force and 
moment in the bending and torsicn equations as 

(2-5a) = La(w,(|!,t) - ^ (M^4)) 


<2-5b) 




m^,(y,t) - m^<w,<|>.t> - 


Incorporation of artificial structural damping g by allowing a complex 


elastic modulus produces the system 
(2-6a) El/l+lg) 0 + " If - % 

(2-6b) GI^(l+lg) - J 1^ + 3^ 3 - a^{y,t> 
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Formal development of the modal equations begins with substitution 
of the Series expansions for w and 4) in terms of the assioned mode 
shape functions and generalized coordinates, With the same number n 
of bending and torsion modes always used, the system order will be 2n . 
Generalized displacements are assumed to have simple harmonic time 
dependence, giving 


n 

(2-7a) w(y,t) = E f (y) q e 
1=1 1 1 


l(Jt 


(2- 7b) 4>(y,C) = E f. (y) q. e 

1=1 "’i 


itJt 


I 

I 

I- 


Generalized displacements for bending modes have units of length, where 
as for torsion they are dimensionless. 

Galcrkln's method Involves substitution of (2-7) into the system 
represented by (2-1), (2-2), (2-3), and (2-4) and then multiplication 
of each term in the bending equation by and each in the torsion 

equation by , followed by integration across the half span. With 

e^^ cancelled the resulting system is 

(2-8a) El (1+ig) E f'”’ f q dy - mu>* E f f q. dy 
* ° 1=1 “i "i • ^ 1=1 "l '^.1 


+ s /**" E f . f q , dy = — F ((|>,w;t)f dy 
e o . - 4). w, ^4), ■' lojt O z ^ w, 

x» I X 1 1 e j 


(1 £ j < n) 
(2-8b) 


GI.(l+ig) r E f!' f. q, 
A 0 4>j 


n 


n 

Z 

i=l 


dy + Jw* E f , f . q, 
4>i <i>^ -t 


dy 


- s /" E f f± q dy = — r~r ® dy 

e 0 4)^ ^iuit o y ' 4)^ 
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The generalized forces can be arranged In terms of dimensionless unit 


generalized forces through the following definitions: 


(2-9a) rv (^,w;t)E dy - 

o z Wj 


(2-9b) 


r " '’'*1 " 1 ii 

Q. 4 -ir ^ Q4 44U.9* ^ 

L 1=1 i -1 ^.J 

r n r. • "I 

. dy = Tlptu^b'ij L Q,, . —7— + S 4 j ^9* I 

[ i=l 1=1 ' J’lj 


list 


The modal integral in the first term of (2-8a), integrated by parts 
twice and with application of cantilever boundary conditions, yields 


/ f"" f dy 

Q »! 


'b 


0 for 


r • dy 
.1 

i * j 


From Appendix A, the i=j term can be expressed in terms of the natural 
frequency of the ^tb assumed bending mode by the substitution 


( 2 - 10 ) 


El r^(£”)^dy = nm>^ /*' 

X o w, w, 0 w, 

j j j 


dy = Jl 
"j 


Similarly, Integrating the first term in (2-8b) by parts and introducing 
Che iiatural frequency to^ of the jth assumed torsion mode leads to 


( 2 - 11 ) 


GIj C fl’ 
d o 4>4 


dy = 


^ f. 

♦3 


(i ^ i) 

(i = J) 


Insertion of (2-9), (2-10), and (2-11) intti (2-8) and f-jrther use 
of modal orthogonality properties gives 
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(1)^ r, 

(2-12a) ^ (1+ig) ^ 


m 

npb- 


Si/ 


e m 


bm Upb*^ 


I'l 


+ E Q , . ^ Q. = 0 

1=1 ^ 1=1 ■1’^'^" '^’l 


(1 1 j 1 n) 
(2-L2b) 


, m _jt ^ _ j m J 

ifpb^ nib"^ ■! U)^ ^({i^ ^ npb^ mb^ ^<J)^ 

s n w n ^w. n 

bm npb- Ij b j+n,i b \i+n,i+n ^())^ = 0 


where the bending equation has been divided through by Hpiii^b^Jl and 
the torsion equation by npu^b**)!. - Inertial coupling, a consequence 
of nonorthogonality between bending and torsion assumed modes, produces 
terms in the modal integrals 

(2-13) i j-^1 t„^(y>f,^(y)dy 


For now it is convenient Co reference tii 

Assumed 


torsion mode natural frequency tii^ 
(Appendix A) are 


(2-lA) 


^,2 = _ 




W. 


An"* 


j r 


Co the fii'St assumed 
mode natural frequencies 


Ratios that are useful are then 
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/ 


(2-15) 


a, 2 

w 

h 


31= 

4 




, ,2 


2 


The transcendental are given in Appendix A for 1 ^ J 5 . 

The equations can he nondimensioralized through the following 
dimensionless parameters, chosen to be as consistent as possible with 
Ref. 1 notation. 


(2-16) 



El 

P = — 5_ 

^ giTF 

a 




Frequency is nondimensionaLized as in Ref. 1 by defining the frequency 
parameter 


(2-17) 


(i+i8)a‘ 5 ti+u)3^^ 


4 ^1 

= IF ^ 

Substitution of (2-15), (2-16), and (2-17) into (2-12) yields the modal 
equations in a suitable form for computation. 


(2-18a) 



(1 £ j < n) 


n 

+ Z 
i=l 




0 
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(2-18h) 




i»l 


"J 


n 

+ E Q 


i-1 J+n,i-Hi ‘’fjij “ ” 


The flutter determinant, of order 2n , results as a necessary condition 
for a nontrivial solution and leads to a complex etRenvalue proVlem to 
determine Z . 

Next, the unit Benerallzed forces Q^.j, must he expressed to allow 
numerical computation. Combining (2-3), (2-4), and (2-7) into (2-9) 
leads CO 

(2-i9a) - npbV <^‘y 


' ^ (y^^A dy> + C It v*bc„a!.-y}2 s f (y)q i 

i=l ^i i 


1=1 j ^i 


(1 1 j 1 «) 


<5,. 


w . n 

X , 


f (y)«iy = e q -rp + e q q ] 

w, i=i 1=1 ^ 


n 

(2-19b) Hpb-u^C / (-M^^ + (is+a)L 1 E f (y)f (y) dy 

\J W W ^ ^ j ^ 

+ - (Js+a)(L +M„) + (J?fa)'L^] E f (y)f (y)q^ dy> 

^ ^ 1=1 1 ^ j 1 

n % 

+ //!ipV^b*a-y)^C Z f" (y)f. (y)-ri dy 

° ‘^1=1 ^ 

n ^w. n 

= npO)*b'*S.t Z 0, . . + Z _ q . ] 


After further reduction, with 
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and carrying out reioalnlng differentiations in y , (2-19) becomes: 


(2-20a) I {L (y)f (y)dy} ~ + (4+a)L 1 • 






(1 < j < n) 


‘ <y)£w^/y>dy + (y)f^^(y)dy]}q^^ 


n n 

^ Q 4 « ^ q* 

i=l ^ i=l <?>£ 


(2-20b) Z 


i=l 


j "j 


J?_ (Z\ 


n 


2Hk> 


+ <^)\l (;>£^ (y)dyJq^, 


n n 

^j+n,i b ^ ^j+n,i+n 


Here, integrals in Che spanwise variable y have been non- 
dimensionallzed. Four different types of modal integrals are encountered, 
including the previously Identified. The three new forms arise 

from the drag coupling terms and are 


( 2 - 21 ) 


fCl) 

ij 

j(2) 

ij 

r<3) 

Ij 


4^i~y)^f;jy)f^^(y)dy 

/^(l-y)fl <y)f (y)dy 

Vi Wj 

j 
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The modal integrals in (2-13) and (2-21) were numerically computed to 
twelve significant digits for 1 ^ ii.1 £ 5 . 

Two additional dimensionless parameters can be introduced, again 
drawn from Ref. 1. 

C 

(2-22) A=5s + 2 C=^ 


The drag parameter C is defined as the applied steady sectional drag 
coefficient divided by the sectional lift curve slope (2f[ for 
incorapresible strip theory). 

After (2-20), (2-21), and (2-22) are combined with orthogonal 
modal integrals recognised, the unit generalized forces are found to be 


(2-23) / 


'j.l 






L (i=l) 

lO (i,^j) 

[L - AL ^ [21. - 4 

w' ij ij 1.1 ij 


[M - AL ]T. . + X 
'■ V w 11 k ij 


-j+n,i+n ~ 


-h M^) -f A Lj (1=J) 

(Wj) 


The depend only upon k, A, and C . Actual computation of the 

Theodorsen function is accomplished by direct use of the ascending 
power series of the modified Bessel functions and , as explained 

in Appendix B. 

Equations (2-4), (2-13), and (2-23) together supply the flutter 
determinant for the VT system using assumed modes, given in Tnatrix 


19 



form In (2-24). Solutions can be found using from one to five assumed 
inodes In each of bending and torsion, and the maximum order of the 
determinant will be 10. 


H[i - n 

-(N n)'^Pi zi + 0 

•A O. n.n 


I _ -MSI, 

I n 'l,n+l in 


' -MS! ,+0 
, ni n,s!+l 


(2-24) 


+ Q, , 

l,2n 

I 

I 

-MSI 

nn 

^^.2n 


Vl.l 




-MSI + Q . , ' 'i(Ml [1- 4-Zl + O . - >0 

nl n+l,n | u 4 n+l,n+l 


-MSI, + 0, , 

In 2n,l 


-MSI + Q_ 
nn 2n , n 


* N 

1 L? Ji. - Cn- )Vl 

I 

i 0 

I 2n,2n 


In seeking solutions, values for M, P, 1^, A, and S must first 
be chosen to specify the wing configuration. Then, for any desired 
value of C and an estimate of reduced frequency k , coisplex eigen- 
values Z can be found from (2-24) via linear matrix eigenvalue techniques. 

A computer program is used to solve for flutter conditions as 
follows. With a first estimate of k chosen large enough so that the 
structural damping, 
lm(Z) 

(2-25) g = j{^(2) 


is negative for all 2n eigenvalues, successively smaller values of k 
are assigned and eigenvalues computed until a positive g is obtained 
for the eigenvalue corresponding to the aeroelastlc mode which encounters 
flutter. Then a zero-finding subroutine locates k for which g=0 
for the flutter node. Dimensionless speed and frequency^ defined by 
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(2-26) 



can easily be computed for the flutter eigenvalue by 
^ “ /ReCZ) 

(2-27) 



The flutter mode shape is also found, since it is related to the eigen- 
vector of the eigenvalue yielding the flutter condition. 

C . Comparison of Modal with Collocation Results 

A rather thorough study was carried out in Ref. 1, covering a wide 
range of practical cowbinatious of the dimensionless parameters. The 
present objective is not to recalculate all of the same data but rather 
to evaluate the effectlverteas of the assumed-mode approach. Consequently, 
(2-24) has been solved at conditions parallel to ones for which results 
are published in Ref. I to offer direct comparison. Also (2-24) is used 
to verify the performance of the VCT system developed in the 
next chapter. 

A comparison of assumed mode calculations with Ref. 1 results Is 

presented in Table. 2.1 for three different configurations. Of these, 

cases (a) and (c) represent a stubby low-aspect-ratio wing (with P = .4) 

whereas case (b) is the opposite extreme of large aspect ratio. A 
El 

typical ratio •= 1.6 , for example, would fix the aspect ratio of 
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cases (a) and (c) at 2 and of case (b) at 20 . The mass ratio para- 
meter of cases (a) , (b) is in the regime encountered by sailplanes 
whereas in (c) it is representative of conventional aircraft. Finally, 
the steady drag parameter is zero in the first two cases but has the 
extrameiy large value C =0.04 in (c) . 

In all cases n=3 already yields adequate modal convergence. 

For the large aspect ratio case the use of one assumed mode (n=l) 
gives a significant error, yet n=2 produces good accuracy. This 
suggests that Che second bending mode is an Important factor in flutter 
of high-aspect-ratio wings, a phenomenon discussed in Chapter 6. 

Additional comparisons between flutter speeds and frequencies 
found by the two methods are offered in Table 2.2, emphasizing their 
relative accuracy as steady drag is increased to the very high value 
C = 0.04 . Evidently good agreement is maintained in the presenc.; of 
drag. 


The mode shape at flutter for two of the preceding cases Is 
presented In Table 2.3 for wings of small and large aspect ratio, each 
with C = 0.04 . Aniplltudes of the ten generalized displacements are 
normalized with respect to |q. | , and the phase angles are referenced 

“Pi 

to the phase of this torsion mode. For P = 0.4 only q, and q 


have appreciable magnitude, whereas the large-aspect-ratio example 

displays a significant contribution by q as well. 

”2 

These flutter mode shapes can be compared to the Ref. 1 results, 
with some effort, as follows. In Ref. 1, the tangent of the angle by 
which the torsional displacement at the wing tip leads bending displace- 
ment in the flutter mode is tabulated; it 1> converted here to an angle 
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In degrees. The spanwise shape of the bending portion of the mode at 
flutter for the two Table 2.3 cases is illustrated in Fig. 5 of Ref. 1. 
Since the flutter inode's bending displacement obtained in that treatment 
was complex, its amplitude at each spanwise station was taken, phase 
differences being neglected in the figure. 

For the present solutions a similar assumption is made to display 
the bending mode shape of the P = 0.004 case in Table 2.3. Contribu- 
tions of q and q have been added vectorically at points along 

Wi 

the span to allow comparison of bending amplitudes with the Ref. 1 
figure. For p = 0.4 , of course, only the first bending mode contributes 
significantly and phase differences are negligible. 

The comparison of bending flutter mode shapes appears in Fig. 2-3, 
where amplitudes are normalized to unit torsional displacement at the 
tip. Phase angles 5 between bending and torsion at the tip are also 
compared. .\c low aspect ratio excellent agreement for phase angles and 
mode shapes is observed, with mild disagreement in the bending mode 
amplitudes. In the P = 0.004 case, for which the Ref. 1 solution was 
made using only five spanwise collocation points, a significant disagree- 
ment in mode shapes and tip phase angles is observed. In spite of this 
discrepancy, the respective flutter speeds differ by only 0.82%. For 
all large-aspect-ratio cases compared this sort of discrepancy in flutter 
mode shapes is observed, where the second bending mode plays a significant 
role and causes appreciable phase differences in bending deflections 
along the span. Mode shapes at flutter are analyzed in more detail 
in Chapter 6. 
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D. Dlveraence Including the Effects of Drag 

One further application of the VT linear model ia predicting 
the effect of drag forces, steady in direction and magnitude and uniform 
across the span, on static aeroelastic stahlUty. Removal of all time- 
dependent terms in (2-1) and insertion of simple strip-theory steady 
incompressible aerodynamic loads gives 
9^W j,2 

(2-28a) El - 2HpV^b(|) + 4 -t (K 4> ) - 0 

X 3y o dy z o 

3 ®(|) 

(2-28b) GI^ + 2npV®b^A<|)^ - = 0 

Subscripts emphasize the.t deflections are static quantities. 

As in the foregoing dynamic analysis, deflections are represented 
by assumed modes and a system of 2n modal equations is derived, which 
has a nontrivial solution only if the determinant of the matrix of 
coefficients is zero which yields the divergence speed with drag effects 
included. For brevity, since the manipulations involved are quite 
straightforward, the final form of the stability determinant is 
presented here. 
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Solution of this determinant for its largest eigenvalue gives the 
divergence speed . This is of interest in the ensuing work, 
particularly for studying dynamic stability of wings having steady- 
state lifting deflections in the presence of drag, which can reduce 
divergence spaed considerably. 

For zero steady drag (C=0) the classical divergence speed comes 
trivially from (2-29), which degenerates for C * 0 . The classical 
divergence mode shape is just the first assumed torsion half-sine mode, 
and the first torsion iriod>l equation uncouples to give the divergence 
speed. 



When A is decreased to zero this classical divergence speed becomes 
infinite; yet when 0 ^ 0 (2-29) will still yield finite solutions 

with A = 0 . 

A direct comparison of solutions of (2-29) with Ref. 1 results for 
the effect of drag on divergence Is presented in Table 2. A. Modal 
convergence is satisfactory, but the modal approach appears to differ 
more significantly from the Ref. 1 analysis for divergence calculations 
than for flutter re.sults. 
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Case 


M =10.0 


a 

S = 0.1 
A = 0.1 
C = 0.0 


P = 0.004 

1 = 0.25 
a 


A = 0.1 
C = 0.0 


i^= 0.25 

S = 0.1 
A = 0.1 
C = 0.04 


(ReE. 1) 


(Ret. 1) 


(Ref. 1) 



% 

2.7175179 

2.7239548 

2.7240004 

1.31052(19 

1.3114559 

1.3114641 

2.7240178 

2.7240199 

1.3114673 

1.3114673 

2.699 

1.309 

4.2621908 

4.0842768 

4.08641(32 

4.0866066 

0.84270? 

0.8849367 

0.8850560 

0.88506.59 

4.0865310 

0.885Q660 

4.032 

0.886 

4.2.60823 

4.260879 

4.260882 

1. 2940 37 

1.2940232 

1.2940250 

4.2608B9 

1.2940236 

4.199 

1.293 


0.239% 

0.00239% 

0.00072% 

0.000077% 


0.918% 


1.337% 


1.4525% 


TABLE 2.1 Comparison of Assumed Mode and Collacation Methods for 
Predicting Flutter Speeds and Frequencies for Taree 
Configurations 
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Case 

c 

Up 


n = 5 

Ref. 1 

% 

Difference 

B 

Ref. 1 

% 

Difference 

P 


0.4 






■M 

mm 

M 

= 

10.0 

0.0 

2. 7240 

2.699 

0.93% 

1.3115 



i 

ct 


0.25 

0.02 

2.7830 

2.753 


1.3071 



s 

= 

0.1 

0.04 

2.S623 

2.823 

1.37% 

1.3024 

1.300 

0.18% 

A 


0.1 








p 

a 

0.004 








M 


40.0 

0.0 

7.1778 

7.127 

0.71% 

0.8651 

0.875 


i 

a 


0.24 

0.02 

6.8588 

6.797 

0.90% 

0.9074 

0.919 


s 

= 

0.1 

0.04 

6.7704 

6.715 

0.82% 

0.9223 

0.934 

1.27% 

A 

ss 










TABLE 2.2 Effect of Drag on Agreement of Modal Analysis with Ref. 1 


Generalized 

p = 

0.4 

P = 

0.004 

Amplitude 

Phase 

Amplitude 

Phase 

Displacement 

'Iw 

1 

0.71276 

219.62” 

1.4515 

195.46” 

2 

0.00132 

21.95” 

0.6288 

217.24” 

% 

0.0000271 

150.13” 

0.00449 

203.35" 

^3 





% 

4 

0.00000408 

147.28” 

0.000674 

-10.94” 

‘’w 

0.00000213 

163.47” 

0.000189 

157.41” 

^5 





’n 

1.0 

0* 

1.0 

0“ 

\ 

0.01503 

222.02” 

0.01454 

147.11” 

\ 

0.000808 

215.78” 

0.01786 

210.39” 

\ 

0.000773 

221.31” 

0.00930 

210.13” 

\ 

0.000186 

217.90“ 

0.00462 

211.73” 


TABLE 2.3 Flutter Mode Shapes for Low- and Hlgh-Aspect-Ratlo Examples. 
(M = 40, 1^ = 0.25, A = 0.1, S = 0.1, C - 0.04) 


27 































n 

C = 

0.02 

C = 0.04 

«D 

% Difference 
From n = 5 


% Difference 
From n = 5 

1 

4.58288 

2.044% 

3.90,105 

1.999% 

2 

4.48660 

-0.0999% 

3.81863 

-0.156% 

3 

4.49174 

0.0144% 

3.82465 

0.0017% 

4 

4.49067/ 

-0.0094% 

3.82403 

-0.0144% 

5 

4.49109 

— 

2.82458 

— 

(Ref. 1) 

4.66 

3.63% 

3.96 

3.42% 



TABLE 2.4 Comparison of Ref. 1 Results With (2-25) Solutions for 
ni'vtargence Speeds in the Presence of Drag. (P « 0.004, 
M = 40., = 0.25, A = O.i). (For C = 0, = 11.1072) 














FIGURE 2-1 
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SECTION A-A 


Physical Origin of 3rag Coupling in Eq. (2-la) 
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Physical Origin of Drag Coupling in Eq. (2-lh) 
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Chapter 111 


development op a general system of equations for 

VERTICAL BENDINC/CHORDWISE BENDING/TORSION (VCT) 

A. Introduction 

Chapter II includes drag loads in transverse bending and torsion 
but allows no chordwlso deflections. In other words, it assumes 
infinitely large bending rigidity in the fore-and-aft direction. In 
this chapter we consider chordwise bending, which is a more complete 
model of the true physical sltu.ation. 

A review of the literature for work cortrernlng VCT motion of 
slender c.antllever beans led to the field of hinge less helicopter rotor 
stability .analysis. The structural modeling of a hlngclcss rotor Is 
essentially the same as desired here, except that the wing has no rota- 
tional velocity. Accordingly, If an adequate model for a hlngeless 
rotor can be found, it can be acLipted for the cantilever wing by rcrwvlng 
the inerti.al and centrifugal tension effects .arising from rotation. 

Much of the work pursued in the hingeless rotor field Includes 
simplifications v;hich cither eliminate or restrict one of the three 
types of deformation, often torsion, in an effort to reduce complexitv. 
This leaves a r<»laClvel? small body of work that treats the full elastic 
problem. A well-known system of linear partial differential equations 
for coupled elastic torsion and bending of twisted nonuniform rotor 
blades is that developed by Houbolt and Brooks (Ref. 6); the iriill-ri 
effort to develop a VCT svsceiri for the cantilever wing centered on 
adapting this formulation. Tn the course of this work, however, the 
elastic coupling terms were found to be insufficient to account for the 
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drag loads In transverse bending and corsion embodied by the terms 
containing In (2-1) „ It then became apparent that this drag 

coupling is actually a nonlinear structural bending- csralsr, effect and 
that the new systeta of equations to be developed should, of necessity, 
retain all nonlinear elastic coupling terms having the same order of 
Importance as these drag coupling terms. 

Further search led to the system of nonlinear cq^lacionB for tvlsced 
nonunifora rotor blades derived by Hodgea and Dowell (Hef. 7). This 
work Involves developmenv of a more complete strain-displacement rela- 
tion chan that of Ref. 6. which Is necessary to obtain the clastic 
bending-torsion coupling terns that produce the desired drag coupling 
effect. The equations in Ref. 7 are valid for straight, slender, 
hoiaogeneous , Isotropic beosei undergoiAfi miMilerate displacittscnts. accurate 
to second order In the sense of a restriction chat squares of bending 
slopes, tvist. and airfoil chord and thickness divided by wing semispan 
are small vich respect to unity. 

Although the final form of the e.jUation« preaented in Ref, 7 might 
appear to be iomediately adaptable to the present cose by setriag the 
rotor rotation frequency to zero and removing Che effects of pretwist. 
this 1 b not entirely true. One important assumption by Hodges and Dowell 
required that the ratio between the transverse and chordvise beam bending 
stiffnesses be a quantity of order one. Whereas this Is a standard 
feature of helicopter blade construction. It usually does not hold for 
conventional aircraft wings of any aspect ratio. As a result, the Ref. 7 
derivation has been carefully retraced for Che specific case of a non- 
rotating cantilever wing having arbitrary bending stiffness ratio. The 
development is outlined in the following section. 
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Development of the Nonlinear Equatrlons of Motio n 


Derivation of the VCT model for the uniform cantilever wing is 
presented here in abbreviated form. Important modifications of the 
Ref. 7 derivation are emphasized, but duplicate manipulations are only 
briefly described. The notation and coordinate system for the canti- 
lever wing is used exclusively. 

The basic ordering scheme presented in equation (4) of Ref. 7 is 
retained. One exception is that spanwlse warping of the cross section 
due to twisting, represented by a warp displacement function which is 
a solution of the Laplace equation over the cross section, is entirely 
neglected here. This assumption is made on the premise that a typical 
aircraft wing airfoil section would have a sufficiently .small thickness 
that warping effects would be negligible within the second-order frame- 
work. The ordering scheme i.<? applied to the energy expressions encoun- 
tered in the variational derivation of the equations, to determine which 
terras should be retained and which discarded. 

The nonlinear strain displacement relations developed In Ref. 7 
have been carefully examined in Che context of the present problem, and 
they are found to apply without modifications. These relations are 
derived from an exact transformation between the deformed and undeformed 
coordinate syteras and originate from the classical definition of strain 
of Novozhilov (Ref. 19) which is based upon increments in the deformed 
coordinates. After approximation to second order consistent with all 
assumptions, the final form in terms of engineering strain and in the 
present notation is 
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f 


u u 

X 2 


Ey ” U* + -y- + ^ u" iCcose + nsino] 


- u”|- ?.slni() + ncos^l 


(3-1) 


Se 


n<f' 


- - f.O’ 


These can be recognized as Identical to equations (24). (2S) , (26), and 
(27) of Ref. 7 after the warp function and pretuisc angle have been 
eliminated. The displacements u^, u^, of the elastic axis and tTc 
principal coordinates s • n of the cross section are illustrated in 
Figure 3.1. 

Develo|Msent of the equations using the indirect method of the 
calculus of variations is based upon Hamilton's principle, which nay 
be stated in the form 

(3-2) / ■ [6(U-T) ■ 6Wldt - 0 

^ % 


The equations are obtained by combining expressions for the first 
variation of strain energy dll , kinetic energy £T , and virtual work 
of external forces <S(!? . 

The first variation of strain energy In appropriate form for the 
standard Euler-Bemoulll beam uniaxial stress-strain relationship is 

(3-3) 6W - ir (o 6 e + o . 6c - + o 6c )df; dn dv 
o yy yy y- >•£; yn 

The first variation of the engineering strains, expressed in terms of 
displacements, is 
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( 

6c “ 6u' + u' 6u' + u' Su* 
yy y X X z z 

- [Ccos(fi + rjalnif)] (Sti^ - u^' 6(()) 


(3-4) < - [^sin(Ji‘+ ricosi))](6u" + u" 5<|)) 

I Z» X 


= n <5 (J)' 

<^e “= ~C S (1)' 

yn 


The 

stresses 

1 are 

! simply 


1 

/ 




1 


a 

= Ee 




yy 

yy 

(3- 

5) < 


r 

y? 

= Ge (- 




o 

+ Ge 




yn 

yn 




Substituting 

/ 


(3-6) < 




(3-1), (3-4), (3-5) into (3-3) yields 

,2 ,2 

n 11 U 

5(1 = ^ JT (E[u' + — ^ - u''(^cosi(> + risinif)) 

area ^ ^ 

- u''(-€sin<(> + ncos<}>) ] [5u' + u'5u' + u' 6u' 

z y X X z z 

- (^cosij) + Tisin<l>) (5u^ - uj^ 6(fi) 

- (-Ssin(}) + ncostfi) (5ii" + u" iS(|)) ] 

z X 

+ G[n^ if 5t^' + (j)’ 6(})']}dn dy 


Rearrangement of (3-6) by grouping of terms having the same virtual 
displacements leads to 
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(3-7) 


6U - / W (5u' + u' 6u + u' 6u‘) + S 54)' 
oyy XX zz y 


+ cosij) + sln4>l(6u” - u” 54>) 


+ [M^ sinij) + cos(};K5u|^ + 64>)}dy 


The stress resultants and moments formed in (3-7), which act on 
the deformed wing as illustrated in Fig. 3.2, are defined as follows 


(3-8) 


V . = JjT dh 


y yy 

^ area 


..2 .,2 


= + “t — + “I ®A^“x 


w y 


where A a. = ff f d£ dr) 
w A 

area 


and iT r, d£! Hr - 0 , by definition 
area 


(3-9) 


= // C a d£ dn 

area ^ ^ 

ut2 ^,2 

° “z ^ ^w®A^“y + 


where ^ .0^ d£ dn 

area 

and ff £ n d£ dri * 0 
area 


^ area 

(3-10) = EI^(u" sin4) + coS(}i) 

where I = ff d£ dh 

X 

area 
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(3-lL) 


Sy = JI (-C O ■*- n 0 -)d5 dri 
area ^ ^ 

= Gi^ r 

where I, = JT + n*)d5 dl) 
a 

area 


At this point in the development the momenta of inertia 

are first introduced. Hereafter, It must be recognized that 1 may 

z 

assume values much larger than . The terms containing the quantity 
e^ , which measures the offset of the tensile axis (area centroid) from 
the elastic axis, will be dropped. The basis for this simplification 

is that, In the final modal formulation of the equations, the retention 

5A e.2 
W A 

of e, terms ultimately leads to a dimensionless parameter, — — , 

rL E<iL 

X 

which appears only as a small quantity added or subtracted with unity. 

Since it will not significantly influence dynamic stability, hereafter 

e, = 0 will be assumed. 

A 

After appropriate integrations by parts within equation (3-7), the 
final form for tne first variation of strain energy, including boundary 
terms, is 


6£J = / {[“(V u')' + f-M cos4> + M sintt>)"]6u 

O yx z 

- (V )'6u + [-(V u')* + sini|) + M cos(^)’')6u 

y y yz z x z 


(3-12) 


+ [-(Sy)' - u^(-M^cos(}) + M^sini})) + u”(M^sin(j) + H^cosi{)) ]6il)}dy 

+ V 5u|^ + >j'V 6ul^ + u’V 5u]^ + S fid'll 
y y'o X y x'o z y z'o y 'o 

i. 5, 

+ [-M cos((i + M sin<i)I<5u'l - [-M cosdi + M slnd)]*fiu | 

Z X )C0 Z X xo 

+ [M sin(fi + M cosd>i'5u'|^ - [M sind> + M cosi|>]'Su 
' z X z'o z X z'o 


37 



Formulation of the first variation of the kinetic energy is greatly 
simplified relative to the Ref. 7 derivation becauae there are no iner- 
tial effects introduced by rotation of the helicopter blade. Since the 
procedure is straightforward and well described in Ref. 7, the details 
of forming the kinetic energy i\c terms of displacement velocities, 
taking the first variation, integrating by parts over time, and expressing 
the resulting form of 6T in terms of time derivatives of , u^, 
u^ , and 4> 3re omitted here. After the ordering scheme has been 
applied, the form of the first variation of kinetic energy, with terms 
retained to second order and corresponding to equation (52) of Ref. 7, is 

6T = jy p^{- [u + $(-£sin4) + ncos(())]5u^ 

area 

- [li^ - $(^cos4) + qsinej)) ] 

(3-2 3) 

- [u^(-^slnij) + ^cos(})) - li^(£cos4) + nsinij)) ] 6<t) 

2 2 

- $[(-£sin4' + ncosij) + (5cos(J) + r)sin(|)) ]6i})}dn dy 

As noted in Ref. 7, the last term in this expression is by defini- 
tion a third' order term, but it is retained In order to include torsional 
inertia in the torsion equation. 

The final form of the first variation of kinetic energy is obtained 
by integration over the sectional area, which leads to the following 
definitions: 
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(3-14) 


1 ■' 
I ‘ 


XT p dC dn 
area 


JT p £ dn 

® area ® 


JT p + n^)d£ dn 

area 


with JT pg n d£ dn = 0 
area 


The final form is 


(3-15) 


6T = r {( -mu + s (tisin(j))6u + (-mii + s $cos<{>)6u 

O X 6 X Z G Z 


+ [s^ sincfi + coscj) - J$]6c(>}dy 


The virtual work of the applied loads is 

(3-16) 6£/ = (F 5u + F 6u + m iS(|’)'ly 

0 X X z z ^ 

Clearly the drag loads will now enter the equations of motion in the 
same manner as the lifting airloads, in contrast with (2-1), since 
and will consist of lift and drag force components. 

Application of Hamilton’s principle using (3-12), (3-15), and 
(3-16), together with (3-8), (3-9), (3-10), and (3-11), results in the 
following quantities being required to vanish: 

Su terms 


u’^ u'^ 

(3- 17a) - [EA^C'i^ + -f- + -f-)]’ 


S T' = 0 
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6u terms 


{-EI^ (~u^cos(j) + u|Jsin<i)) cos(f 


(3-l7b) + El (u"alu<li 4- u"cos(j>)sln(t!}” + mu 

X X ^ z X 


- s^(|) sin<J) = 0 


6u tem^s 
z 


{eI (“u" cos 4> + u” sin.(t!)ain(i) 
Z X z 


(3-17c) 


+ EI^(u^ sln{|> 4- coscf))co3(|>}" 4- mu^ 


s $ cos<|) - F =0 
e z 


5(1) terms 


(3-17d) 


(GI, ({)')' - u"[- EI (- u*' cos(|) 4- u" sin(|))cos4) 
M z z z 


+ EI^(u^ sinij) 4" coz<ti)siTii))] 


4- u” [EI. (-u” cos(|» 4- u" slnt})) sini}> 4- EI (u” slivl) 4- u" cos$)cosil>] 
X Z X Z X X z 


- £ li sln<t) - s ii costfi + J«ii - M “0 
ex e z ^ ^ y 


In (3-17a'), since It is known that for the nonrotating contilever 
wing the spanwise tension T will he zero everywhere, the expression 
can be Integrated leaving 
u'^ 

(3-L8) u*4--^4-~-=0 
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(true only with - 0) . This constraint indicates that the only 
manner by which axial deflections enter the problem is by the 

purely geometric dependence on u^ and due to shortening induced 

by the lateral deflections. 

In replacing sln4) and cos<}) by small angle approximations, it 
would be consistent with the ordering scheme to approximate the cosine 
by unity. Due to the possibility of a large ratio of bending stiffnesses, 
which can result in u^ being large compared with u^ , it was found 
necessary to keep the second-order approximation of the cosine at this 
step. Tnis is done to derive properly certain elastic coupling terms 
involving and , while maintaining symmetry in the inertial 

and stiffness matrices of the final matrix equations. 

Substitution of 


(3-19) 


cosf 1 - 
sin4) “ (|) 


into (3-17) leads to 


(3-20a) [- EI^(- u"(l-((i^) + u^(J>) + EI^(u^<{»^ -t- u"(ti)]‘’ 


+ mu - s od) = F 
X e X 


(3-20b) 


[El (- u» + u"(!))c|. + El (u"0 + u" (!-(!)='))]" 

b X Z XX z 


+ mil' - s $ = F 
z e z 
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- (Gl.f)' - u"[- + u"((.) + El (u"(|)=' + 


X X" 


(3-20c) 


+ + uy^) + El (u"(|) + u"(l-(|.^))] 

A Z X Z X X z 


-su(()-su +J(fi = jn„ 
ex e z y 


Reorganlzaclon of terms and Introduction of the displacements w 
and V for and u^ , respectively, gives 


(3-21a) 


(3-21b) 


<3--21c) 


{■[EI^il)^ + EI^(l-it)^)]w" - (EI^ - EI^)v"(t)}'' 


+mw-s5=F 
e z 


{- (EI^ - EI^)w"(|) + [EI^(l-(j>^) + <|)^ET^]v’'}" 


+ mv-S(f(|> = F 
€ X 


Gl^cj." - (EI^ - EI^)[(w"^ - v"^)(}) - v'V (1-^)1 


+ svd) + sw-J<£+ tn„ = 0 
e e ^ y 


The underlined terms are reasoned to represent higher-order effects and 
are dropped. In the case of the inertial terms coupling chordwise 
bending and torsion through the offset s^ , this is a third-order 
effect relative to terms like mv . The final form of the three-DOF 
system of equations, with some regrouping, is 


(3-22a) 


{El w" - (El, - El )(v"<|) - w"(|)^)}" 

X Z X 


4* raw - s !|i = F 
e z 
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(3-22b) 


{EI^v" ~ CKT^ - El^)(w"(ti + + mv = 

(3-22c) GI^i|)" - (EI^ - EI^)[(w"^ - ^)(J) - v'V] 

+ sw = J$+m = 0 
e ^ y 

Thts system is elastically and inertially self adjoint, which 
assures that the stiffness and Inertia matrices eventually assembled 
during modal analysis will be symmetric. The tarms containing 
(EI^ - EI^) represent the nonlinear coupling between the torsion and 
bending degrees of freedom; all remaining stiffness, inertial, and 
applied-load terms represent the same familiar forms encountered in 
linear beam theory. Terms in (3-22) which are underlined do not appear 
in the Ref. 7 equations and are retained here as a result of the absence 
of a restriction on the bending stiffness ratio EI^/EI^ . Strictly 
speaking, when this ratio is large compared with unity, the single 
underlined terms will increase in relative importance while the double 
underlined terms are negligible in magnitude. 

The nonlinear equations of notion (3-22) are next adapted to 
permit analysis of stability about a steady-state deflected position 
due to an equilibrium lifting flight condition, which could be level 
flight or a steady pullout at a high load factbr. Small time-dependent 
perturbations about the equilibrium operating condition are Introduced 
by expressing the deflections w, v, and ifi in terms of steady-state 
equilibrium deflections ‘J’o small unsteady perturba- 
tion quantities '?X* 
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(3-23) 


' w(t) w + w, (t) 
o 1 

( v(c) = V + (t) 

O i 

^ <|)(t) = <))^ + 


First, Che steady equilibrium deflections only are substituted 
into (3-22) to obtain a nonlinear system of equations for the equilib- 
cium solution. The resulting nonlinear steady system is 

(3- 24a) {ei w" - (El - El )(v"i!i - = F 

XO 2 X OO 0^0 2 

0 

(3-24b) {EI v" - (EI - EI ) (w’’c|) + = F 

zo z X o o 0 0 X 

o 

<3-24c) GI,4i" - (EI - EI )[(w"^ - - v" w"] + m =0 

' ' d^o 2 X O 0 ^0 oo y 

■' o 

Appropriate Steady lifting aerodynamic loads are inserted for L and 

o 

M, , and the assumed steady drag, which entered equations (2-1) in an 
entirely different manner, is Introduced through L 

Next a linearised system of equations in the time dependent small 
perturbation deflections is obtained by substituting (3-23) into (3-22), 
subtracting the nonlinear equilibrium equations, and discarding products 
of the perturbation quantities. The linearized unsteady equations of 
motion are 


(3-25a) {El w" - (EI, - EI ) [i(i v” + (fi.v" - - 2(ji w"(j!,]}" 

X 1 2 X O 1 1 O O 1 


+ iri)^ - 3 - F =0 

1 el 
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(3-2 5b) 


{EI v'' - (El - El )t(J) wV + w"^ + + 2<p v"4ij}" 

z 1 z X ‘^0 1 1 o 0^1 ^ 


+ mv, - F “0 
1 x^ 


(3-25c) 


Gid^*’ - (Ei^ - Ky {[2v;w-- - 2v%*']4>^ + [«;;2 - 


- w’’vV - v'V.'} + SW, -J$,+m =0 

0 1 o 1 el 


Analysis of stability with this system can be done with the standard 
techniques for linear systems, but first the coefficients must be found 
for a given equilibrium flight condition by solving (3-24). The loads 
appearing in this equation, including the chordwise forces, must be 
expressed as linear functions of perturbation displacements and 

(j)j^ . No dependence of aerodynamic loads on the fore-and-aft motion 
Vj^ will be considered in this analysis. 


C. Comparison of the VT and VCT Models 

The VCT model (3-24) and (3-25), although different In appearance 
from the VT model (2-1) , in fact reduces Co the same form when 
T ->■ oo and the steady lift is zero. To Illustrate this, first imagine 
that the steady loads applied to (3-24) are 
f 

F - 0 
z 

o 

(3-26) < m =0 

] ^ 
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where D is the same assumed drag force, constant in magnitude and 
direction along the span, as is considered In (2-2). The solution to 
(3-24) immediately gives w^ = <})^ •» 0 , leaving just the chordwise 
bending equation as 

(3-27) El v’”' = D 
z o 

Substitution of xEI for El in (3-25) together with 
X z 

w = (j) =0 gives 

0 0 " 


(3-28a) 

(3-28b) 

(3-28c) 


El w'"' - (T-l)EI (c|),v")" + mw, 
X 1 X ^1 o 1 



= 0 


tEI vV 

X 1 


+ mv. 




(t-1)E1 (-V 
x*^ o 


.. 2 , 


I, - v"w"] + s W- - 4- m 

1 o 1 el 1 y 


= 0 


The equation for uncouples and can be disregarded. Integration 
of (3-27) twice, using the zero shear and moment boundary conditions 
at the free end, gives 


(3-29) 


tEI v" = D 

X o ' 

= - M 

z 


where (2-2) has been used. Substitution of (3-29) into equations 
(3-28a, c) yields 

(3-30a) El wV” + (M «,)‘' + mif - s (J - F >0 

3C JL i Z X .L 6 X Z ^ 
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(3- 30b) 


T-1 


GI .d)'' - M (wV + v"4.J + s « - + TO +0 

d^l T z 1 o 1 el y. 


Finally, requiring that T , with the consequent vanishing of , 
causes (3-30) to reduce exactly to (2-1), 

For the TOOte general case of a steady lifting condition charac- 
terized by 0 and 1 * 0 , (3-24) and (3-25) can still be cast 

into a fom similar to (2-1) , in which the displacements become 

dependent upon and since x ->■ » . First, the nonlinear steady 


equations 

(3-24) are rewritten with replacing 

El 

2 

together 

with some 

rearrangement of terms to give 



(3-31a) 

El w"" - TEI ~ [v" - w"4) )<t ]" = F 
xo XT o ooo z 

0 



(3- 31b) 

TEI [v" - ^ (w"4> + v"(j)^)]'* = F 

X O T O O 0 0 X 

0 



(3-31c) 

GI A" - TEI - v")w" - ] + 

d^o T X 0^0 0 o 0 

m 

V 
' o 

= 0 

Again, as 

T ^ the elastic bending curvature about 

the 

n principal 

axis (Fig. 

3.1) of the airfoil section should go to zero. 

But now, the 


section is displaced to a position fixed by the deflections w 


and <P of its elastic axis, and the true elastic curvature about the 
o 

n axis is now recognized as v” - w”$ 

In equation (3- 31b), which is the fore-and-aft equilibrium equation, 

as T becomes large and v becomes small, the tex'W v"6^ is of 
^ o o o 

higher order and can be neglected relative to , Removing this term 


and letting 


T-1 


1 leaves 
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(3-32> 


X o 


o o 


F 

X 

o 


Integrating this equation twice, using the shear and raonient boimdarv 


conditions on v and w at the free end, leads Co 
p o 


(3-33) tEI [v" - w’'0 1 = F f - M 

X O O O X z 

o o 


Here it is assumed for convenience that F is constant along the 

X 

span; of coarse the integrations could as well be p'erformed for any 

Known spanwise variation of F . M is similar to M of (3-29), 

' o o 

except that it represents the bending- tors ion coupling effect upon the 

moment-curvature relation. As t — ^ the ouantirv (v" - w"(J ) must 

go to zero according to (3-33). In the limit the v deflection becomes 

o 

dependent upon and as a result of bending- tors ion coupling. 

The two remaining Indepes’.dent equations in (3-31) are already 
arranged so that the coupled curvature quantity in (3-33) can be 
recognized. Sutstitution for this quantity leads to 


(3-34) 


Elw“" 

o 



-^ol* 

o 



Cl.r' - — IM w" 
d^o T z o 
o 


V 4* 1 + m “ 0 
i> O V 

■ o 


As T , with v*’({) recognized as a hiisher-order term, the aonl inear 
o o 

steady equilibrium equations finallv become 

(3-35a) EIw"" + (M 1" - F 0 
o z o z 

o o 
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(3- 35b) 


GI .lb" - M w" + m 

d z o y 

o o 


0 


These resemble (2-1) wlrh dynamic terms removed. 

A tinal exercise is to demonstrate that the linearized unsteady 
perturbation equations (3-25) also reduce in a systematic fashion to a 
VT system when When (3-25) are rearranged with tEI^_ 

replacing EI^ and the combination of steady deflections in (3-33) 
replaced by , the linearized system becomes 


T-1 


(3- 36a) El w"" + [M <!),)" - tEI ~ [v'.' - 4) w'.' - w"(l),)it) )’ 
X 1 ‘ Z ^1-' X T 1 o 1 o 1 


+ nw, - s Ifi, - F “=0 
1 1 z^ 


(3- 36b) 


IKI [v” - w” + + (J)2v" + 2 ^ v>,)V 

X 1 T ^O 1 0^1 1 ^o o’l 


+ mv, - F = C 
1 X, 


(3-36c) 


GI.<t)'’ + tEI — [(vV - 4.,w" - vj"4i,)w'’ + 2v"iJ) v/' + ' 

d 1 X T 1 11 o 1 o o o 1 o 1 


M w + s w - j 4> 4- m =0 

z 1 G 1 ^1 y, 

O . -^1 


The underlined terms are of higher order for t > 1 and are discarded, 

T-1 - 

When T is sufficiently large that — ~ 1 » equation (3-36b) becomes 

C3-37) 

i 
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The quantity B « defined as 


(3-38) B = v!,’ - (j) w',' - , 

1 o 1 o 1 


is recognized as the time-dependent curvature, in the linearized 
perturbation deflections, about the r| axis of the steady-state deformed 
cross section. Thus it in expected that as t ->■ “ the curvature B 
must approach zero. 

The dynamic linear moment M is defined as 
(3-39) ^ - tEI^ B 


Equations (3-37) and (3- 36a, c) can 'be restated, after use of (3-38) 
and (3-39) together with t , as 


(3-40a) El wV + M (fr, + [M + mw, - s - F = 0 

X 1 z 1 z, o 1 e l z, 

o 1 1 


(3-40b) M" = F - mv, 
"l ’'l ^ 


(3-AOc) - M w” - M w*’ + s w - -f m =0 

d I z L z, o el 1 y.. 

0 1 1 


Since B 0 . the acceleration term in the second equation could be 
expressed entirely In terms of w and <{) accelerations by working 
from integration of 


(3-41) 


V” - 


fo'^1 


w"$ 

o 


1 


in this linearized unsteady system, M (y,t) represents the moment at 
station y due to the instantaneous chordwlse inertial loads and applied 
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linearized unsteady drag forces acting outboard of the station. Because 

of the dependence of M upon w, and <t>. , equation (3-40b) cannot 

1 i 

be uncoupled to allow for a separate calculation of M" . Thus, even 
though can be eliminated, (3-40) still involves three unknowns 

"l , ‘f’l > \ ■ 

The purpose of this section has been to demonstrate analytically 
the connection between VCT representation of the cantilever wing, given 
by (3-24) and (3-25), and the VT system given in (2-1) and used in 
Ref. 1. It can be concluded that the VCT model remains valid for 
arbitrarily large bending stiffness ratio T , and that the drag 
coupling effect has been saCisfactorally accounted for by the nonlinear 
elastic bending-torsion coupling terms. Since the forms obtained in 
(3-35) and (3-40) are but special cases of (3-24) and (3-25), actual 
solutions will be found, using the latter system only, for practical 
values of T . 
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Chapter IV 


FLUTTER VELOCTTY CONSIDERING VERTICAL 
BENDING/ CHORDWISE BENDING/TORSION (VCT) 

A. Modal Equations for Small Oaclllatory Motions About a Steady-State 
Deformatio n 

The model for torsion, transverse bendint;, and chordwlse bending 
of the uniform CL'-ittlever wing developed In the previous chapter is 
used to determine flutter velocity uith of the same sec of assumed 
inodes used in (2-1). Use of ti assuraed inodes In torsion, n In trans- 
verse bending, and n in chordwlse bending results in a set of 3n 
modal equations in terms of modal generalized displacements. i>lnce the 
assumed modes satisfy the natural as well as geometrical boundary 
conditions, which were obtained during the application of Hamilton's 
principle, Galerkin's method can be employed to transform the equations 
into algebraic relations In the generalized displacements. The nonlinear 
steady-state equilibrium equations (3-2A) become nonlinear algebraic 
equations in the steady-state generalized displacements, which are 
solved iteratively. These displacements determine the coefficients of 
the linearized unsteady model by applying Galerkin's method to (3-25). 
Then the velocity is determined for which simple harmonic motion of 
this system is possible (neutra] dynamic stability). 

The steady aerodynamic loads for the equilibrium equations <3-24) 
are specified in terms of incompressible strip theory. A typical air- 
foil section (Fig. 4-1) has its zero-lift line inclined to the free- 
stream velocity V by the angle a + . The resultant steady 

lift L (acting in a direction perpendicular to V ) and the moment 

d. 

o 

m are given by 
0 
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(A-1) 


L = 2HpV^b(a + <}) ) 
o 

< 

m = 2lTpV^b^ (-.+ 3 ) (a + ) 

o 


A drag force D is also present. These forces must be resolved into 
the axis system x » y, k fixed with respect to the wing rone. The 
required transformation is 


(4- 


( F =1 cosa 4- Dsma 
I % a 
I o o 

< 

F = Dcosa - L sina 
X a 

o o 


Assuming sina ~ a and cosa ~ 1 and neglecting the z component of 
drag shows that the steady aerodynamic forces to be used with (3-23) are 


F = 2Il0V^b(a + (^ ) 
z o 

o 


(4-3) / 


K = 2IIpv^b(C - - a4> ) 

X 

o 


I ra = 2npV^b^A(a + iji ) 

V ' o 


The drag force represented by C is constant spanwise, and the 
definitions in (2-22) have been adapted. 

The steady-state deflections are now expressed in terras of the 
assumed modes (cf. Appendix A for definitions) 
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( 4 - 4 ) 





The indicated differentiations in the first two equations are carried 
out, and appropriate Integrations by parts of the resulting modal 
integrals is performed. The first terms of each equation are then 
expressed in terms of modal natural frequencies as in (2-10) and (2-11) 
which leads to the form 


, tq*’ - (T-1) E S II , q* q" 

( j J U=1 V=1 ■*' % li 


El 


n n 


(4-6a) ( 


El 


n n n 


F- ^ ^ijUM 

U=1 V=1 1=1 ^ i 




- 2npV^bK Z I, q? = 2npV=‘b«l 
, <p 

V=1 •' 


2B a 


El 


n n 


Tm U)" I q” - -rr (T-1) E Z H „ . q' q, 


(4-6b) 


w -v f 

j j 


, , VUJ (p w 

v=l !i=l u 


El n n n 

- -rr (T-1) Z Z E R., q" q” q“ 

' U=1 P=1 i = l \ ''i 

r 2B. 

+ znpv^bto Z I q“ = 2noV*blL - (c - a*)I 

v=i J % L .1 


(1 < j < n) 


El 


n n n 


T nl + -7^ (T-1) Z Z Z R, 

4.J 2 <{.. t 


^ ^ ^ Vil 'C '^w 

v=l h=L i=l ^ ^ V U 


ET 


n n n 


(4-6c) 


-yr (T- 1 ) I I Z R ,, q* q' q* 
V=1 U=1 1=1 ^1 


= 2RpV^b^e.A 


1 

n(j-'s) 


and are properties of the assumed bending modes described in 

Appendix A, and represent the modal integrals previously encountered 

and defined In (2-13). The nonlinear bending-torsion coupling terms 
give rise to two new forms of sodal integrals. 
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(4-7) 


«juv = (5)dy 

J IJ ^ 


^ijyv = ^0 ^w 

-' 1 J ^u ^v 


Equations (4-6) can be arranged in the final format used for computa- 
tions by dividing the torsion equation by HpV^Z-b^ and the two bending 
equations by UpV^jlb , using modal natural frequency relations in (2-14) 
and (2-15). and nondimensionalizing with the parameters given in (2-16) 
and (2-26) . 


MP-! ^w. (T-l)MPi 

4V,- w 1 ' ' a 


n*N 


(4-9a) 


\ "wsv"'* b 


u 


j U" b U‘ 

h»i V=1 i=l \ V=1 


4R. 

^ ^ k ( “21 

^ ' i 


tMPI ^v. 
/ n^N! — ^ -r^ 
I J D 


KPl 

- (T-1) 


(4-9b) 


n n “w 

I E H ^q'’ 

, , vpj A' b 

v=l y=l ■' V 


1 


n n n v. ( n 

+ 2^2 R-..,. q! 2a s l, q" - tr* (c-«^> 

, , , , ijpv ^<1 b , Iv^d) N,n 

v=l y=l 1=1 -* v=l j 


4B. 


(i j j« n) 
/ 


V'inHi-hr -77T- - A]q« 


. MPi f n n n m? ^w 


O O 

q q 

n n n 

y y y U V 

b b "juu b b 

v=l U=1 i=l ■' i u=l u=l 


n n 

z I -a. 


\ V, 


2A a 


ncj-y 
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A solution to (4-9) can be found for a given wing configuration 
once drag C , root angle of attack a , and dimensionless speed U 
have been specified. The solution procedure involves Newton's method 
and is described in Appendix C. 

When simple harmonic motion Is assumed, the linearized perturbation 
e>|uatlons of motion (3-25) are converted Into linear algebraic equations 
via the same steps described in Chapter II In deriving nodal equations 
(2-18) from (2-1). The perturbation displacements are expanded in 
terms of n assumed modes and generalized displacements, 

w (y.t) - 5: f (y)q^ 

i-l 1 1 

n 

(4-10) V (v,t> »* I f (v)q c 

i-i “i ''i 

<>^(y,t> - 1' (y)q^ 


where and q have dimensions of length, whereas q. arc 

dimensionless. Anticipating solution for flutter bound.iries using the 
V>g method, structural dampinp. g is included by introducing a complex 
clastic modulus, and the staple hamonic airloads .are assumed to have 
been expressed la terns of unit generalized forces as In (2-9). With 
the addition of the fore-and-aft bending degree of freedom, the 
generalized forces now appear as 


(4-Ua) 
(14 J in) 


n n 

- r n , ♦ .> O, |e 

i-l *» 1-1 i--^2n C, 


Iwt 
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(4- I lb) 


/ F (C,.w ;c)f dy 
o 1 1 


n n 


- Hcw‘b’t( Z Q. 
1-1 


+ Z o. 
i-1 


If 

'1 


iw't 


(1 ^ J ^ n) 

(4-1 Ic) m (0..W ;t)f dv 

o 1 1 


n 

• n£x-j*b'«i z 

i-1 


*'l+2n.l 



n 

T 0 

1-1 


*1. 

'1 


}c 


U;i 


The ooission of colusms In the attav h^ivinf; n +l ^ I ^ 2n 

) » * 

reflects the jssuaptlon th.it there Is no dependence of unste.idv .ilr- 
loads upon fore-and-.if t ewtlons . 

After C.ilerkln's aethod h.is been laple«i>nted for eqixitlons 
by the suae aanlpul.it ions re<;ulred to develop (2-17). the svsten of 
3n aodnl equations which determine linear st.ibilitv about the stendv- 
equilibriua deflections appear in the fores 


f {M - n- N 


a 


(4-12a) 




n n 


^ ^ Nll.v < ■=; > IT 

1-1 •' U-! v—1 ’0 V 


1-1 li-l 


1 


n n 

+ r. U^-i)MPi, 7-1 - H 
l-l 


*^v It n 

U 


zi :: H, , - 2 r : r ^ ,q* 

vi • iui b , . 1 7. b 

u-l li^l V— 1 U 


I 




S9 



n n ’v 

j i-l ® U -1 ^ 


{i-l2b) 




+ {M - Tri*'M‘!MPl Z) ^ 
j a b 


n n 


+ r {(T-1)H?1 Z l I qt q! }-— 

, , 3 , , tiuv i> ft b 

i“l Vi"l V *1 " U V 


n n 


* /, ‘<t-nw‘„ n r h + 2 z z . -^1 

U-I U"1 V«1 ■’ 11 


* ‘*5+n.i+2n^**$j " ^ 


(i < ] *•' n) 


( '■ 


n n 

I {(:-l)HPi^ 7.\ I H, 
I 


,, n n ’w 

“-22 2 R q" -^] 


a • -, -uil b ~ "uivi ■’ft. b 

u-l u»l v»l • V 


i -SmT * Q ^ } — ~ 


n n 


V V. 


<z-. 2 c> ^ * r .(T-n«pi,zi r -jE * 2 r r 

j 1»I U-i •’ u»l V«1 V 


q q ft o 

n n V n n ^ 

■* r, {Ct-DMPI ZI r ^ R ,.-r^-rr- z Z R _ u y, 

u-l v»l U=1 v=l b b 


n 


l-I 


* - "’U-y-zlh* ■ 0 


Here widal Integrals defined in (4-7) 

and (4-S), which appear In this systeo oniltiplled by steady-state 
generated displacements. The cof ficiencsi of terms which couple the 
steady deflections into this systent must first be constructed by the 
indicit^T^rd siaanations be fere the eigenvalues can be determined. 
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The unit generalized forces required in (4-12) will now be expressed 
in terms of the two-dimensional incompressible unsteady aerodynamics 
adopted for (2-2) and (2-3). It will be necessary, as it was for the 
steady aerodynamics, to transfer the airloads from wind oriented axes 
to the coordinate system fixed at the wing root. For the present case, 
the circulatory part of the unsteady lift will be assumed perpendicular 
to the direction of the freestream velocity, whereas the noncirculatory 
portion will be assumed to act normal to the chord of the airfoil sec- 
tion in its steady-state deflected position. The Inclusion of unsteady 
leading-edge suction effects will be considered in the next chapter. 

To separate the circulatory and noncirculatory contributions to 
the lift, in (2-2) and (2-3) can be rewritten in the form (cf. Eq. 

(4-126) of Ref. 5) 

w. i4>-, 

(4-13) La(w^,<!)^;t) = npb*w2[-^ + k 

+ npb^u^C(!c) ^ V ^ ^ X 

Here the first term on the right represents noncirculatory lift L 
and the second term the circulatory portion L . Referring to 
Fig. 4-2, these terms can be expressed as resultant forces in the xz 
frame by 

( F = L cosa h h cosd) 
r = - L sina + L sin(|) 
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KakinK the assumptions cosa ~ 1 and cos<^^ ~ 1 yields tha Identity 
(4-15) r- - L (v,,0,;t) 

7 a X X 


when Ic is assumed Chat sina ~ a and slnd ~ 4> > the chordwise 

o o 

force components become 

( 4 - 16 ) = npbVa C(k)l-^^ - ^ ~Y 


w 

+ nob^w*4> t-r^ + 

o b k 1 


itot 


Preparatory to steps that lie ahead, it can be seen that when the 
unit generalized foro7s are determined from (4-11) with the modal series 
for inserted, modal Integrals of Che forms 

4 

1. j V 

4 

i 3 V 


will be encountered as a consequence of the noncirculatory contribution 

ho (4-16). Since these integrals do not occur elsewhere and the effect 

they Introduce is expected to be minor, Che noncirculatory contribution 

to r is neglected altogether. This simplification has the result 

Chat the unsteady noncirculatory force illustrated In Fig. 4-2 acts In 

Che z direction, at all spanwise stations. 

The remaining circulatory contribution to F can be expressed, 

*1 

using the notation in (2-3) and after some manipulation, as 
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w 

(4-17) = npb’u^al(l-Ly) - ^ + ^ + A(1-L^)1 ^j 

3> 


With Che aerodynamic loads given in (4-17), (4-15), and (2-2), 
where t) = m^(Wj^,<J)j^;t) , tlie unit gecterallzed forces can 

be c jvelopedc starting from (4-11), by means of the same steps followed 
In assembling (2-22). The final result is 


r ‘>J,i ■ ('•v 


(i=j) 

(i#j) 


^1,i+2n = - (L(^ - 


'fj+n 


T 


(1-L,)a (,i=j) 

(ij^j) 


(4-13) 

(l<j<n) 

(Ki<n) 


-VI h.i 


'>1+2i>,1+2«-|>5[«^ - + M„) + a\) 


0 

'V,i+n 


= 0 for 1 < \* < 3n 


(i“1) 
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B. Solution Procedure for FluCCer Velocity of a Llftin;; WitiR 


The 3n medal equations (4-12), with (4-18) inserted, can be 
expressed in matrix form 

(4-19) ([Mgl + [Q]){q} - ZtkgKq} = 0 

The mass and stiffness matrices H and K are real and svmmetric 

s s 

and are written out in Fie. 4-3. The terms In K include sums of 

s 

products of modal integrals and the steady equilibrium generalized 
displacements found as solutions of (4-9). The aerodynamic matrix 
Q , whose elements appear in (4-18), Is complex and nonsymmetric . 
Equation (4-19), therefore, represents a complex eigenvalue problem 
for the complex frequency parameter Z . Its solution yields damping, 
speed, and frequency as in the case of (2-25) and (2-27). 

The logic used to compute neutral stability conditions from 
(4-19) is diagrammed in Fig. 4-4. The primary difficulty encountered 
when steady deflections are introduced is that a preliminary estimate 
of speed must be made before the eigenvalue problem can be solved. 

Steady deflections for are used to generate coefficients in , 

and unsteady aerodynamic loads are then computed for large enough 
reduced frequency that the eigenvalue corresponding to the aeroelastlc 
mode which flutters has negative (stable) damping g . Successively 
smaller values of k are then substituted and aerodynamic terms 
recomputed. The eigenvalues are recalculated until positive damping 
occurs for the flutter mode. From this point a zero-locating sub- 
routine determines the reduced frequency for which damping g is zero, 
whereupon one can calculate flutter speed . The entire procedure 
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is repeated with a new estimate until the speeds and Up 

are adequately matched. The flutter speed determined for zero lift 
serves as a Rood first estimate for U^ . With care, close aRreeraent 
of the two speeds can be achieved in three or four iterations. 

Three assumed modes in each degree of freedom, corresponding to 

a system order of 9, are found to give adequate convergence for all 

cases. The model integrals were computed to allow n < 4 , a task. 

which required numerical integration of 100 quantities of type 

and 40 of type H. . 

ipv 

C. The Nonlinear Elastic Coupl ing ier nis 

The need to retain all terms in (3-22) , including third degree 
nonlinear, in order to model adequately the nonlinear elastic hending- 
torsion coupling mechanism is now demonstrated by means of typical 
applications. Of course, neglecting higher-order nonlinear effects 
would have the appeal of reducing complexity. For example, removal 
of all third-degree nonlinear terms would result in the elimination 
of all terms containing quadruple model integrals from both 

the nonlinear steady equiflbrium system and the linearized dynamic 
stability analysis. In order to examine the effects of such approxi- 
mations, numerical experiments were conducted wherein higher degree 
terms in both the steady and unsteady modal equations were neglected. 

First, static deflections are considered. Clearly, If the 
equations do not adequately represent the steady-state deflections 
of the wing over a realistic range of lifting conditions, then any 
linearized unsteady perturbation analysis based on the same level of 
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approxlnuttlon cannot be expected to succeed. The nonlinear steady 
equations (3-23) eij.l be solved at three levels of approximation: 

(1) Linear terms only retained. 

(2) Second degree monllnear terms Included. 

(3) All nonlinear terms included. 

The linear set of modal equations is just (4-9) with the nonlinear 
terms removed: 


(4-20) 


rN. 


MPi 

a 



n 


2 I 
v=l 





n^N*! T 


HPi 


V . 

b 


+ 2a 


n 

E I. q* 

v=l 


4B. 


nn. 


(C - a*) 


r->n^(j-5.0^ 




2Aa 

WTj^) 


The torsion equations are now seen to be n uncoupled and iramedlately- 
solveable relations for the n generalized displacements q® . The 
results permit solution for the bending displacements, thus eliminating 
the need for matrix operations. Since all elastic coupling terms are 
absent, deflections and ((»^ are independent of the stiffness 

ratio T and the steady drag parameter C . A very important conse- 
quence is that the mechanism by which drag influences divergence is 
missing. One. concludes that (4-20) is a particularly poor representation 
for .steady deflections in the presence of steady drag. 

Retention of the second-degree nonlinear effects from (4-9) involves 
adding to (4-20) all terms containing the triple modal Integrals 
H. . The system becomes fully coupled, and solutions are found by 
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the sasue iterative scheme described in Appendix C for solvioR the full 
nonlinear system. 

The behavior of the second degree nonlinear solution is compared 
with Che linear deflection from (<i-20) in Figs. 4-5 and 4-6 . Here 
it is shovm how the deflection in the first vertical bending mode and 
the first torsion mode vary with Increasing speed for a fixed angle of 
attack a = 0,01 radlons of the wing root. The parameters K , P , 
i^ „ A , and T for these examples correspond to the idealization of 
a sailplane wing discussed in Chapter VI, and sufficient modal con- 
vergence Is assured by using n=3 . Numbers on the ordinate of Fig. 4-5 
show the actual twist in radians of the wingtlp due to the first torsion 

mode; the vertical deflection in semichords of the wingtlp due to q“ 

''l 

Is just twice Che value read from the ordinate of Fig. 4-6. 

The important influence of drag on the steady deflections Is 
evident. But there is poor correspondence between the second degree 
nonlinear deflections and the divergence speeds indicated in Fig. 4-5 , 
which are solutions of Che linear static stability determinant (2-28). 
This discrepancy is most pronounced when C=0, and the rapid divergence 
of the second-degree nonlinear deflections for all values of drag 
suggests trouble with the second-degree approximation. 

Figures 4-7 through 4-10 show solutions of the complete non- 
linear syetem (4-9) for the same wing conflguretion, together with the 
linear results. Correlation with the divergence predictions appears 
to be excellent, and the sudden blowing up of deflections characteristic 
of the second-degree nonlinear solutions is not encountered. Since the 
iterative procedure used to solve for nonlinear deflections commences 
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with the linear solution as an initial estimate, nonlinear solutions 
which do exist for C=0 above the classical divergence speed cannot 
be obtained. 

The full nonlinear solutions conform with nonlinear behavior 
expected by intuition. For zero drag, the vertical bending and torsion 
deflections should fall below the linear solution owing to the effec- 
tive increase in stiffness "seen" by each degree of freedom due to 
defcrination in the other. This effect is observed. The slight rear- 
ward chordwise displacement (Fig. 4-9) for C=0 when and 

are large comes from bending-torsion coupling, and the negative contri- 
bution of the second torsion !Sode (Fig. 4-10) reflects a redistribution 

of elastic twist toward the wingtip where the curvature due to w is 

o 

less. Drag alters the deformed state radically and causes large 
displacements at much lower speeds. 

A clue to the reason why the second-degree terms alone are inade- 
quate is found by looking at the sensitivity of solutions to changes 
in the bending stiffness ratio T . The second degree nonlinear 
solution for bending deflections, given in Fig. 4-11, reveals increas- 
ingly poor behavior as T is increased. Conversely, the complete non- 
linear solution behaves as intuition would anticipate, becoming 
insensitive to changing T as this parameter grows toward “ (Fig. 4-12) . 

In Chapter III the behavior tor T-^“ of the steady nonlinear 
equilibrium equations was analytically investigated. In the discussion 
following equation (3-33) it was pointed out that the quantity 

represents the actual curvature about the principal axis 
of chordwise bending for the. deflected airfoil section, and that this 
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quantity should vanish as t . In this llinit any chordwlse 
deflection Is a geometric consequence of elastic coupling between i}) 

o 

and ; terms which contain the product approach a 

finite value. Now, referring to equations (3-24), it can be seen chat 
In the second degree approximation this product is retained in the 
second equation but not In the vertical bending equation (where 
is dropped but '''q'Pq retained). Likewise In Che torsion equation, 
w"** is dropped but v"w" is kept. The result is that terms 
containing rv^ remain and blow up in the limit t . Accordingly, 
the conclusion is reached that, for structures representative of air- 
craft wings for which T is reasonably large, the third-degree terms 
containing the product must be retained. It may be added that, 

although this elastic coupling effect was neglected in the helicopter 
blade equations developed in Ref. 7, in that case the ordering scheme 
required T to be on the order of unity. Another point Is chat the 
remaining third-degree terras In (3-24) , which contain the product 
, actually can be neglected when T is large. 

Since it has been determined chat third degree terms cannot be 
excluded from the steady nonlinear equilibrium analysis, it follows 
that they must also be kept in the linearized unsteady pf-rturbatlon 
equations (3-25). The original dynamic stability equations used early 
in this investigation, on the other hand, consisted of a linearized 


unsteady perturbation system based on the second-degree nonlinear 
approximation. That is, they are the result of removing all terras 
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as was an Extreme sensitivity to large T which could not be justified 
physically. Furthermore, linear steady-state deflections were used — 
a particularly unsuitable approximation for Inclusion of steady drag. 

Examples of these original flutter calculations appear in Figs. 

4-13 and 4-14. The former gives a comparlscn of the early results 
with the full analysis by (4-12) tor the sailplane example, showing 
how the simpler approximation differs significantly in flutter speed 
even in the presence of moderate steady deflections. Figure 4-14, 
based on a different wing configuration, gives an idea of the difficulty 
encountered for large T when using the simpler analysis. (Note chat 
the steady deflections are quite small in view of the fact that this 
exaisple Is a large-aspecc-ratio wing). All of these stability boundaries 
abruptly terminate, at which point the eigenvalue solutions began to 
behave erratically. In contrast, flutter solutions of (4-19) can be 
obtained for arbitrarily large steady deformations — indeed for tip 
deflections well beyond practical limits of material linearity. 

D. Comparison of VT and VCT Calculatlotis 

The flutter velocities found using the VCT model can be checked 
directly against velocity computed from the VT model developed in 
Chapter II. Although such a comparison restricts the former method to 
the special nonlifting case = 0 , Its basic approach of 

linearized perturbaticn motions about a steady-state deflection solution 
can nevertheless be tested by inclusion of steady drag. This is 
because the drag effect enters the stability determinant through 
coefficients depending on . The cbordwise bending equation of 
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notion is uncoupled when v ■ v * 0 . and fluccer nodes involve 
just vertical bending and torsion. 

The VCT scabillty equations were analytically reduced to the VT 
form (2-1) in Che last chapter by specifying n ■ 0 and t -* •» . This 
agreement is reflected in numerical results, as evidenced by Fig. 4-15. 
For an except ionnlly 1,-irge magnitude of drag. C » 0.04 . chosen to 
magnify the importance of drag coupling and hence the steady deflections 
In the stability analysis, the two methods are used to compute flutter 
spf^eds for the same nonliftlng conf Iguration. The VCT analyses were 
mpde for x ranging from 1 to 10,000. 

As is always true for EI^ = EI^ , the drag coupling effect 
vanishes at T = 1 with flutter speed unaffected by drag. At Che 
Other extreme, for x = 10,000 the computed flutter speeds for 
C = 0.04 differ -by a mere 0.052%. The Ref. 1 flutter speeds calculated 
for this same configuration are also shown and agree quite well. 

When one examines the VCT dependence of flutter speed on T 
(Fig. 4-15), the effect of drag is apparently insensitive to T over 
the range of this parameter representative of practical aircraft 
applications. Hence the VT system (2-1) does not suffer by its 
inherent assumption that X . 

Table 4.1 presents a further comparison of results from the two 
approaches for a different configuration, in which modal convergence 
Is emphasized. Of course, since X = 50 exact agreement cannot be 
expected . 

Convergence of flutter speeds and frequencies and also mode shapes 
for Che VCT solutions is documented in Table 4.2 for zero steady 
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lift and In Table A. 3 for a lifting condition specified br ci • .01 
rad. In the latter case, the apparent slower convetitcnce is related 
to the steady lifting deflections (also tabulated), vhldit also v&ry 
with n . 



n 

Up. W 

i Bp. vcr 


2 

4.26M23 

, 



4.25S457 


3 



4.250351 


r 

4.2M)$a2 


4.25813S 


5 

4.260SS9 

L ■ 


4.258316 

_ 


TA«tC 4.) Coatparlson of Mo«ial Fluttor Sp««<i« Conpured with the 
VT «nd VCT 5yt«t«tts. <M - 40. . P - 0.4. I - 0.25. 
A » O.l. S - O.i. T - 50. C - 0.04) “ 



n ■ 1 


n ■ 2 


n ■ 3 


n ■ 4 



Up -4. 

15027 

Up - 4.183899 

Up - 4. 

183883 

Up -4. 

183916 


“ 0. 

85254 

flp « 0.0S876R 

Op » 0. 

88757 

Hp ■ 0. 

8875S 

■ 

MIPL. 

PHASE 

AMPL. 

PHASE 

AMPL. 

PHASE 

AMPL. 

PHASE 



224.54* 

2.0650 

223.87* 

2.0659 

223.87* 

2.0659 

223,87* 

1 










- 


0.6345 

-46.85’ 

0.6309 

-46.84* 

0.6308 

-46.85’ 










'*w 

- 




0.0107 

- 4.46* 

0.0107 

- 4,26’ 

\ 

- 


- 


- 


0.0022 

1.89° 

\ 

1.0 

0* 

1.0 

0’ 

1.0 

0’ 

1.0 

0* 

\ 



0.03067 

19’'. 53“ 

0.03075 

199.70“ 

0.03075 

199.70’ 

\ 

- 


- 


0.00340 

- 5.55* 

0.00339 

- 5.67’ 

\ 

- - 


- 

i 

- 


0.0059 

220.49* 


TABLE 4.2 Modal Convergence for Flutf.er at Zero Steady Lift. (M = 9.4, 
P = 0.01, = 0.25, A =» 0.1, S = 0.1, T = 25., C = 0) 
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n - 2 

n = 3 

Op - 3.5495 

Up = 3.5329 

Op - 0.68200 

Op = 0.67468 

1.19C55 

1.18142 

C. 01336 

0.01333 

- 

0.000817 

0.002270 

0.002296 

-0.000279 

-0.000284 

- 

-0.000033 

0.0087100 

0.0087591 

0.0000373 

0.0000915 

- 

0.0000017 

1 "* 
AHPL. PHASE® 

AMPL. 1 PHASE® 

2.4492 I 227.45 

2.4765 1 227.68 

0.52487 j -26.65 

0.30948 j -26.27 

I 

0.01637 j - 6.54 

1 

1 

0.84595 j - 4.77 

1 

0.85914 j - 4.89 

0.03071 j 179.03 

0.02979 1 180.02 

' 1 

0.00923 j 172.67 

1 

1 

1,0 j 0 

1.0 1 0 

0.07178 -2X,65 

0.07798 1 -20.83 

- 

1 

0.02793 - 4.83 
1 


n = 4 

Uj. = 3.53X5 
= 0.67276 

1.18074 

0.01333 

0.000818 

0.000151 

0.002305 

-0.000284 

-0.000033 

-0.000005 

0.0087656 

0.0000983 

0.0000090 

0.0000013 


AMPL. 

PHASE® 

2.4860 

227.76 

0.51068 

-26.23 

0.01612 

- 6.57 

0.00250 

- 1.47 

0.86569 

- 4.98 

0.02916 

180.20 

0.00532 

172.76 

0.00122 

172.56 

L.O 

0 

0.07947 

-20,73 

0.02925 

- 5.06 

0.00998 

- 8.60 


TABLE 4,3 Modal Convergence for Flutter at Steady Lift for a = O.Ol rad. 

(M = 9.4, P = 0.01, i = 0.25, A = 0.1, S = 0.1, T = 25 , C = 0) 
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FIGURE 4-3(a) Mass and Stiffness Matrices Appearing in Eqnation 
(4-19) 
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FIGURE 4- 3(b) 




FIGURE 4-4 Flow Chart for Solution Procec 
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Using V-g Kethod 

















FIGURE 4-6 Amplitude of the First Vertical Bending Mode. (Same Conditions as in Figure 4-5) 
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FIGURE A-8 Amplitude of the First Vertical Bending Mode, (Same Conditions as in Fig. 4-7) 
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Figure 4~10 Amplitude of Che Second Torsion Mode. (Same conditions 
as in Fig. 4~7) 
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85 


of the First Vertical Bending Mode vs. T, Using Second-Degree Nonlinear 
lar Approxirnatlons . (Remaining Parameters Same as In Fig. 4-5) 






nonlinear 
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FIGURE 4-12 Amplitude of the First Vertical Bending Mode vs. T, Using Full Nonlinear and Linear 
Approximations. (Remaining Parameters Same as in Fig. 4-7) 
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FIGURE 4-13 Flutter Speed as a Function of Steady Wing-Tip Displacement for Two Different 
Approximations in the Analysis. (M= 9.4, F = 0.01, 1 = 0.25, A = 0.1, S = I 
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Chapter V 


THE DETERMINATION OF AEROELASTIC MODES FOR ARBITRARY VELOCITY 

A. Incompressible Strip Theory Alrloada for Arbitrary Motion 

The procedure for determining flutter velocity developed in 

Chapter IV has several drawbacks. It requires a matched point analysis 

in which an estimated velocity used to calculate steady deflections 

has to be iteratively matched to the lowest calculated flutter speed 

U_ for the proper structural damping. Consequently the solution for 
r 

a stability boundary over a range of lifting conditions can be lengthy. 
Furthermore, intermediate computations have no physical significance 
and are jf limited qualitative value. The behavior and degree of 
stability of individual aeroelastlc modes, which becomes more inter- 
esting with the addition of the fore-and-aft bending degree of freedom, 
has proven to be difficult to deduce from the V-g solutions. The 
only quantitative information available pertains to the neutral stability 
conditions found for the mode which experiences flutter. 

We would like to have physically meaningful information regarding 
dynamics of the system at any desired speed. That is, we would like 
to know the complex eigenvalues of the aeroelastlc modes at subcritical 
and supercritical velocities. Obviously a major drawback of the V-g 
method is its dependence upon simple harmonic air loads. 

Several solution procedures were studied, which replace the 
Fourier transformation with respect to t by Laplace transformation 
with respect to t and obtain at least an approximation to the modal 
stability below and above the flutter velocity. The p-k method, a 
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British flutter analysis technique (Ref. 8), uses the same simple 
harmonic airloads but assumes that these loads are approximately com- 
plex eigenvalues s = a + iu , where Ui is the frequency used in 
calculating the airloads and |oj « |w| . 

Another approach, commonly applied in helicopter blade stability 
analysis, involves use of quasisteady aerodynamic theory by assuming k 
is small enough to permit C(k) ~ 1 . Laplace transformation of the 
system then yields only linear and quadratic terms in s . Linear 
matrix techniques can then be used in determining the roots. This level 
of approximation neglects entirely the effect of the unsteady wake upon 
the circulatory airloads and is not appropriate for the magnitudes 
of reduced frequencies observed in many flutter calculations by the 
V-g method. 

A third possible course is to apply an augmented-state method, 
which approximates the actual unsteady aerodynamic loads for arbitrary 
motion with a transfer function relating airfoil displacements to loads 
having a rational Laplace transform, resulting in a linear matrix eigen- 
value problem for the aeroelastic modes. Goland and Luke (Ref. 9) 
used this route to study wing bending-torsion flutter. They adopted 
the R.T. Jones (Ref. 10) approximation to the Wagner indicial lift 
function to express unsteady airloads in rational form, taking the 
Laplace transform in time. In addition to their accurate description 
of the basic bending-torsion aeroelastic behavior at all flight speeds, 
Goland and Luke demonstrated that the severity of flutter cannot be 
reliably Inferred from solutions by the V-g method. 
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Each of the techniques mentioned above attempts to gain information 
about aeroelastic inodes whose elgeiiveluea have nonzero real parts. The 
true effect of the unsteady wake on the aerodynamic Loads for arbitrary 


motions is approximated to varying degrees. This Is done quite well 
for most motions in the case of augmented-state methods, marginally 
in the p-k method, and not at all in the quasisteady case. 

For present purposes all of these schemes were rejected in favor 

* 

of the more exact approach developed by Edwards (Ref. 11 ) . An impor- 
tant contribution of Ref. 11 (adapted from Scars, Ref. 21j is the 
definition of a generalized Theodorsen function to represent the exact 
circulatory two-dimensional incompressible unsteady airloads in the 
Laplace domain for arbitrary motions. The generalized Theodorsen func- 
tion is expressible in terms of the modified Bessel functions of complex 
argument and as 


(5-1) C(i) 


Kj (s) 


K (s) + K.(s) 

O 1 


where 


sb 
S = V 


Although previous investigators had recognized that this form was 
convergent for the right half plane, representing divergent oscillatory 
motions with Re(s) > 0 , the restriction on the integral definitions 
of the modified Bessel functions caused some investigators to believe 
that convergent oscillatory motions (s in the left half plane) could 

*See also Milne (Ref. 22). 




ii Miiiiii; 


iiak' . 
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not be So represented. Edwards observed that K (s) and K, (s) are 

O 1 

defined and analytic throughout the s-plane except for a bri.Tich point 
at the origin. Wien one places a branch cut along the negative real 
axis, C(s) can be shown by analytic continuation to relate circulatory 
loads and displacements throughout the s-plane except along this cut. 

It is, in effect, an "aerodynamic transfer function" in the Laplace 
domain . 

With substitution of s = ik in (5-1), the fam5.1iar Theodorsen 
function of reduced frequency for simple harmonic motion is recovered. 
Although arbitrary motion is now being considered rather than simple 
harmonic motion, the two approaches yield similar forms when the initial 
conditions arising in the transfotms are neglected. In fact, the simple 
harmonic airloads (4-18) and (2-3) can be used for arbitrary unsteady 
motion simply substituting C(s) for C(k) and s for iuj . 

The modified Bessel functions are computed from their ascending 
power series expansions, as mentioned in Ref. 11 and described in 
Appendix B, Since the transforms of aerodynamic loads will be multiple- 
valued functions because of the branch point of C(s) at the origin, 
the convention 


- n < Arg(s) _< H 

is used for the cut on the negative real axis. The generalized 

i0 

Theodorsen function is computed by Edwards in the form C(fe'^ ) for 
representative values 0 _< 6 ^ II . It is shown to approach 1 as 
r -»• 0 and h as r -»■ <» for all 6 . 
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with the ability to compute unsteady two-dimensional incompressible 
airloads for arbitrary motions in hand, the simple harmonic stability 
analysis developed in Chapter IV can be generalized for this case. All 
of the true aeroelastic eigenvalues and eigenfunctions can be obtained 
for any prescribed speed . Stability can be displayed with root 

locus diagrams. The solution technique is developed in the following 
section. 


B. Solution for the Aeroelastic Roots b y Means of Assumed Modes 

Formally, the procedure for developing the modal equations needed 
for the true aeroelastic modes throughout the complex plane begins with 
Laplace transformation of the linearized equations of motion (3-25) . 

The transformed perturbation displacements are then expressed as series 
expansions in the assumed modes as 

n 

w (y;s) = I f (y)q <s) 

1“1 1 1 

(5-2) V (y;s) = I f (y)q (s) 
i=l ^i ^i 

^ (y;s) = I f. (y)q, (s) 

^ 1=1 ^1 ^i 


A system of homogeneous, linear, algebraic equations in the generalized 

displacements q , q * then be derived by Galerkin's 

'^i ''i '^i 

method, as before- Nontrivial solutions are given by the zeros of the 
determinant In s . Since the coef ficieuts in this determinant which 
arise from the aerodynamic loads will contain the nonratlonal function 
C(s), this will not be a polynomial eigenvalue problem. Roots 
s = a + ioj will thus have to be obtained by iteration. 
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Since full development of the modal equations is analogous to that 
for the simple harmonic case, the stability determinant in a will 
simply be constructed directly from the modal equations for simple 
harmonic motion defined by (4-12) , (4-18) , and (2-3) . One replaces lu 
by s , ik by s , and eliminates g . The elements of the aerodynamic 
matrix given in (4-18) for simple harmonic loads become, for arbitrary 
motion. 



(i=j) 

(i/i) 


in which 


(5-4) 


L = 1 + 2 
w s 

L, = ^ [1 + ZC(i)] + 2 

5 S 
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These aerodynamic loads are based on the assumptions discussed prior 
to Eq. (i-i3). 

With the same mass and stiffness matrices M and K illustrated 

s s 

in FiR. 4-3, and the aerodynamic matrix Q whose elements are defined 
by (5-3), the matrix form of Che lolal equations in s becomes 

(5-5) lUM^l + (ql) IK^]}iq(s)} = 0 


(5-5) may he compared with the simple harmonic form (4-19). It is 
convenient to define a dimensionless Laplace transform variable 


(5-6) 



which is related to the reduced Laplace transform variable through Che 
dimensionless velocity by 


(5-7) 


sb _ o 
V ' U 


The stability determinant thus takes the form 
(5-8) Ip^lMj + [Q]) + iK^ll = 0 


Zeros of this determinant will yield 3n exact roots for the aero- 
elastlc modes in terms of the 3u assumed modes. These roots describe 
modal frequencies and stability at the speed used to calculate the 

steady-state deflections which enter as coefficients in the stiffness 
matrix. 

A computei: program was developed to locate numerically the zeros 
of the determinant (5-8); the logic is outlined in Fig. 5-1. This 
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format proved to be quite convenient for constructing root locus 
diagrams having either the speed U , the root angle of attack a , 
or the drag parameter C as the changing parameter. Any of the Impor- 
tant aeroelastlc modes, at any degree of stability, could be traced 
through the complex plane as long as the ilni>.ial guess of s was 
sufficiently close to Its particular root locus. 

Since the determinant order is 3n and n > 2 is desirable to 
model adequately the physical system, numerical expansion of the 
determinant was not practical. A library subroutine, employed to 
calculate the complex determinant, proved to be the source of numerical 
difficulties. It was found that for 1 ^ a ^ ^ and a = 0 the 
Fig. (5-1) program gave accurate results when compared to parallel 
V-g method neutral stability computations, agreeing to at least seven 
digits In flutter speeds. For n * 2 and a 0 , which gives rise to 
steady deflections due to lift, similar good agreement was encountered. 
But for n ^ 3 and a ^ 0 , the program converged on zeros which 
did not match the neutrally stable V-g predictions and were obviously 
incorrect from a physical standpoint. Subsequent investigation revealed 
that the numerical difficulties originated in the library subroutine. 

When a = 0 the fore-and-aft bending degree of freedom is dynami- 
cally uncoupled, and the order of the determinant which was actually 
computed by the library subroutine was reduced to 2n . Thus, for 
n = 3 and a = 0 , the actual computed determinant was of order 6, 
whereas for n = 3 and a 0 the order was 9. For the latter case 
the actual magnitude of computed determinants was often 0(10^^) , 
while for a = 0 the n = 3 determinants 0(10^*^) . For n «■ 2 
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and a = 0 determ;i.nants again were 0(10^®). The magnitude of the 
computed determin^inc thus appears to be related to the numerical 
difficulties. Since n - 2 results are judged to model the problem 
adequately and never encountered numerical problems, correction of the 
above difficulties was not pursued. 

As a result, all root loci shown herein for steady lifting condi- 
tions involve n =• 2 . As will be discussed, however, this restriction 
doe.s not compromise Che modeling of the physical system nor prevent 
qualitative understanding of its behavior. Furthermore, numerical 
results always agree reasonably well when compared with V-g computa- 
tions for a = 3 , a 3 ^ 0 . 

The algorithm used to estimate the zero.s cf the determinant in 

the S plane, using an initial guess , is illustrated graphically 

in Fig. 5-2. The complex determinant is first calculated at and 

at the two related points s + .001 and s + 1.001. Points A and C 

o o 

are then determliied, at which linear extrapolation in the two orthogonal 
directions predicts that the real part of the determinant v/ill vanish. 
Similarly pcints B and D are predicted, for which the imaginary part 
vanishes by extrapolation. A new guess for the root s^ is then 
determined as the intersection point of dashed llne.s in the figure- 
The process is repeated until satisfactory convergence ia realized. 

This simple scheme worked quite well and never failed to converge on a 
root, usually within four or five iterations. The convergence criterion 
used was generally |ls^| - |s^|| £ 10 ^ . Typical performance of the 
algorithm is documented in Table 5.1. 
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C. Mode Shapes for Aeroelastlc Modes 


The mode shapes associated with roots determined for arbitrary 
motion by the determinant iteration method were conveniently calculated 
with the same linear eigenvalue routine applied in solving the simple 
harmonic stability problem. For a root computed by iteration, the 
nonrational aerodynamic tsmts containing C(s) in the matrix 0 
can be immediately evaluated. Thus one Is led to a conventional matrix 
eigenvalue problem, which ootains the same root as one of its eigen- 
values and also provides Its eigenvector. This approach worked well 
and, as a bonus, verified the accuracy of the roots computed by the 
determinant iter.ation scheme. 

Since numerical difficulties were encountered with the determinant 
evaluation routine for n = 3 and a 0 , this eigenvalue approach of 

recheckinp its results for n - 2 , a 0 (and for all n with a = 0) 

Is valuable. It offers the only means of verifying computed roots 
lying off of the iui axis. Correlation to at least six significant 
digits was always observed. Moreover, the accuracy of the n - 2 . 
a ^ 0 determinant iteration solutions has been checked for a few 
representative cases by letting n •= 3 in the eigenvalue routine, with 
the known n = 2 root as a first guess to evaluate the aerodynamic 
loads, and iterating until the true n - 3 root is obtained. Two 
cases for which this was done, together with their eigenvectors, are 
shown in Table 5.1, where n = 4 roots are also given. The n = 2 
results obtainable by the determinant iteration routine are thus seen 
to be acceptably accurate even for the higher frequency mode. 
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0 . Inclusion of Unsteady Chordwise Loads Due to Leading Ed^e Suctlo ii 


Since motion of the ulng in fore-and-aft bendlne is permitted, 
unsteady chordwise loads can participate In the dynamic stability 
problem. Unsteady two-dimensional Incompressible airloads given In 
Eqs. <A-18> and in (5-3) are strictly based upon the assumption that 
the Instantaneous resultant unsteady lift on anv airfoil section along 
the span is always perpendicular to the direction of the free stream 
velocity. Two-dimensional incompressible potential flow theory, however, 
does predict an unsteady leading edge suction force which arises from 
Che Inverse squat 2 root singularity of the vorticity distribution along 
the airfoil chord at its leading edge. This effect will be included 
into the analysis within the. framework of Che linearized unsteady per- 
turbation theory used to determine stability. The effect of the unsteady 
propulsive force on stability can then easily be isolated by comparison 
of roots computed for the (4-18) airloads directly with roots determined 
with the airloads derived in this section. 

The existence of a leading-edge suction force due to the leading- 
edge sirigulatlcy was determined by Yon Karman and Sears (Ref. 12). 
Greenberg (Ref. 13) in developing the propulsive force on an airfoil 
in an oscillating stream, scares chat a propulsive force acting on the 
airfoil in the upstream direction. 

(5-9) F = HoC^ 
s F 


arises from the unsteady vorticity distribution which behaves at the 
leading edge as 


(5-10) 


Yl 


-1 


2C, 


/b 


-1 
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In this application the coefficient will contain contributions 

from the superposed vorticity distributions due to Che steady and Che 
urstendy lifting flow fields. Therefore Cp will be a sum of steady 
and unsteady parts, and C“ will involve a steady portion, a linear 
unsteady cross multiplication part, and a nonlinear quadratic unsteady 
terra. In the context of the linearized unsteady perturbation approach 
to stability analysis, only the cross multiplication term will enter 
the dynamic equations. To include consistently the nonlinear unsteady 
propulsive force effect on stability, the nonlinear structural coupling 
terms discarded during linearization would have to be reintroduced. 

As a consequence of linearization, the unsteady propulsive force 
can be included only when both steady and unsteady vorticity distribu- 
tions are present. Thus the case of zero steady lift will have no 
contribution due to this effect to the state of stability. The effect 
will become increasingly pronounced as Che steady lift is increased. 

The vorticity singularity strength in (5-9) was given by 

r 

Garrick (Ref. 14) for an airfoil oscillating in a uniform stream as 
(5-11) Cp = {[h + Va + ba(4-a) ]C(k) - Jjbct} 

with h positive downward. Converting to present notation, introducing 
the generalized Theodorsen function, and introducing superposition of 
steady and unsteady deflections gives 



(5-12) 


Cj, = {(a-ht^)V + [-w^ + + b(Js-a)ij]C(i) - 

The propuislire force , using (5-9), is 
(5-13) Fg = 2npb{(a+<{.^)^ + 2V(a-Ki^) [ (-w^ + 

+ b(*5-a)<j)^)C(s) - 'ib<|)j^] + [(-Wj^ - tfijV + b(l 2 -a)ij>^)C(s) 

- 

The last squared term Is the nonlinear time-dependent contribution and 
is neglected hereafter. The first squared terra is the propulsive force 
on a flat-plate airfoil at incidence in steady flow. The resultant 
steady aerodynamic force should act at right angles to the free stream 
velocity in potential flow, and this steady propulsive force can be 
interpreted physically as the component which tilts the resultant lift 
vector, obtained by summing the pressure distribution at right angles to 
the chord, forward to become normal to the alrstream. 

The propulsive force can thus be seen to correct for the chordwise 
component of the lift which is computed normal to the airfoil chord. 

The assumption incorporated into (4-18) that the unsteady circulatory 
lift acts at right angles to the airstream must be discarded and t.he 
force assumed rather to be normal to the chord of the airfoil in its 
steady-state orientation. The assumed direction of forces is shown in 
Fig. 5-3. 

In the structure axis system, the loads are 
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(5-U) 


p T. _ 1? cosd) + L Sind) 
s o fl o 


F “ L cosi}) + F aini 
.'1 o 3 o 


The followlns? assumptions are made 

cost}) ~ 1. 
o 

sin<t» ~ d> 
o 

F ainib << L 
a a 


so that 


F " L 

a 


F - - F + L il) 
Xj s a o 


The lift for arbitrary motion in time ta given hy Ref. (5) as 
(5-lb) = 2nobV[- Wj, C(s) + V C(s) + b(?i-n)4)j^ C(s)] 


+ - ba^^l 


Combining (5-16) and the linear unsteady part of (5-13) into the second 


of (5-15) yields 


(5-17) = -4npVb(cx + 'i(()^)(- Wj^ C(s) + <J>j^ V C(s) 


+ ((-*5-a)C(s) - ^)b«i>j^] + ITpb®V(})^(t»j^ 
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The loads can now be arranged into the unit generalized force 
format of (5-3) by Laplace transformation in time, substitution of the 
assumed modes as in (5-2) , and formation of the generalized forces from 
relations such as (4-11). After this work, which is straightforward, 
the elements of the aerodynamic matrix which incorporate the linearized 
unsteady propulsive force are 


(5-18) 




+ (2-L ) E H. 


• ) ^ n ,, ql 


(l<i<n) 
(IS j<n) 


= 2[L, - - 4 + A(l-L )]I,^ a 

j+n,i+2n (1> s w ^ ji 


+ [L, - 1 - -f + (2-L )A] Z Y. ^ q° 
^ ® v=m+l 


where the new modal integrals appear. 


(5-19) Y.^^, = / f f , f . dy 

ivi 0 w, q) <P . 

■J 1 '■V 


The remainder of the terms of the aerodynamic matrix remain the same 
as in (5-3) . 


Prior to actual calculations a further approximation is made. 

The terms in (5-18) which depend on sums of q? are neglected, 
eliminating the need to compute the • This is equivalent to 

assuming that the lift L in Fig. 5-3 is aligned with the z-axis and 
that is neglected relative to a in the linear term of (5-13). 

Strictly, this simplification will alter the results somewhat, hut It 
is not expected to change the overall effect of the propulsive force 
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on stability and does simplify computations. The first order trend of 
the effect of 100% leading edge auction on stability is the main point 
of interest and should not be affected. 

The changes made to the aerodynamic matrix are therefore substitu- 
tion of the terms 


(5-20) 


Q 


j-Hn,i 


*^j!-n, i+2n 


I 


2(1-L )a 
w 


0 


(i-j) 

(1'l‘j) 



- h - i + A(l-L )]!,, a 
s w ji 


(1< j<n) 


(1 < i < Ki) 


for their counterparts in (5-3) . 

The program described in Section B includes the option of using 

either of these unsteady aerodynamic force systems, avid a comparison of 

their relative effect on stability is made in Che next chapter. Except 
2 

for the - — term in the second of (5-20), incidentally, the newer 
s 

system simply involves doubling the magnitude of the terras in (5-3) 
that are replaced. 
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M = 40, P = 0.005, = 0.25, A = 0.1, S = 0.1, T = 25, a = 0, n = 3, 

C = 0 , U = 6.5 


ITERATIOK 

s 

DETERMINANT AT s 

0 

-0.072968 + 10.0582071 

0.326x10^ + 10.1266x10^'^ 

1 

-0.078290 + 10.067052 

-0.259x10^ + iO. 504x10^ 

2 

-0,079541 + 10.065839 

0.106x10® + 10.504x10^ 

3 

-0.079527 + 10.065910 

0.106x10® - 10.178x10® 

4 

1 

-0.079526 + 10.065910 

-0.281x10''^ + 10.499x10^ 


M = 40, P = 0.02, = 0.25, A = 0.1, S =-0.1, T = 25, a = 0.02, n = 2, 

C = 0, 0 = 7 


iteration 

s 

DETERMINANT AT s 

0 

-0.0055179 + 10. .119612 

0.261x10^^ - 10.253x10^° 

1 

-0,0014167 + 10.1074866 

0.125x10^° + 10.428x10^® 

2 

0.0005501 + 10.1085285 

-0.459x10^ + 10.135x10® 

3 

0.0005027 + 10.1084861 

-0.129x10® + 10.178x10® 

4 

0.0005026 + 10.1084861 



TABLE 5.1 Performance of the Determinant Iteration Algorithm for Two 
Cases, One Konlifting With n = 3 and One at Steady Lift 
With a = 0.02, n = 2 
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STEADY-STATE DEFLECTIONS 


M - 40 

P - 0.005 

i = 0.25 
a 

A = 0.1 
S = 0.1 
T » 60 

a * 0.01 RAD. 
C - 0 
U = 7 


ROOT FOUND BY DETERMINANT 

ITERATION WITH n = 2: p = 0.03579 + lO. 55875 


<C 

V 

q; /b 

<‘/b 

q" /b 
0 ^ 

^ o 

<^2 
_ o ^ 


L_ 


n = 2 
1.97768 
0.022566 


0,003170 

-0,00363 


0.0068192 

- 0.0001112 


2.01496 

0.023221 

0.001453 

0.003534 

-0.000403 

-0.000044 

0,0072034 

0.0000129 

-0.0000322 


2.019088 

0.023274 

0.001464 

0.000271 

0.003594 

-0.000407 

-0.000046 

-0.000007 

0.0072456 

0.0000303 

-0.0000130 

-0.0000103 


RESULTS OF LINEAR EIGENVALUE ANALYSIS FOR n = 2, 3, 4: 



n = 2: 
p = .035795 

+ i. 55875 

n = 3: 

p = .036828 + 1.52903 

n = 4: 
p = .036486 

+ i. 52309 


AMPLITUDE 

PHASE 

AMPLITUDE 

j PHASE 

AMPLITUDE 

PHASE 


4.2245 

1.2064 

215.40'“ 

-40.51“ 

4.6573 

1.1808 

1 217.75“ 

' -37.53“ 

[ 

4.7497 

1.1815 

218.17“ 

-37.07“ 


- 


0.03432 

-18.25“ 

0.03544 

-17.97“ 


- 


- 

1 

1 

0.00530 

-12.46" 


0.083839 

- 7.71“ 

0.87304 

j - 8.79“ 

0.88887 

- 9.12" 


0.05128 

181.23“ 

0.05146 

181.49“ 

0.050487 

181.93’ 

</b 

- 


0.009128 

1 169.96 

0.009314 

169.12" 

i 3/b 

- 


- 


0. 00 2 135 

165.92" 

- 4 

-*2 

-*3 

1.0 

0.05780 



0.0“ 

-33.85* 

1.0 

0.07738 

0.02511 

I 0.0” 

-30.30“ 

1 -19.79“ 

1 

J 

1.0 

0.08213 

0,02969 

0.00912 

0.0“ 

-29.93" 

-18.96" 

-22.80" 

L_ 


TABLE 5.2 Modal Convergence for Convergent and Divergent Oscillatory 
Aeroelastic Modes (see Fig. 6-24(a)) 
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ROOT FOUND BY DETERMINANT 

ITERATION WITH n = 2: p = -0.216505 + 11.632293 

RESULTS OF LINEAR EIGENVALUE ANALYSIS FOR n = 2, 3, 4 



n = 2: 


n “ 3: 


n = 4 : 


I 


p = -.2165053 + 

il. 6322931 

p = -.21485 + 

11.59221 

p = -.21564 


11.57794 


AMPLITUDE 1 

PHASE 

AMPLITUDE ' 

PHASE 

AMPLITUDE 

r 

1 

PHASE 

or 

0.02478 ' 
1 

175.84“ 

0.02627 [ 

175.18“ 

0.025894 

I 

I 

174.95" 

q /b 

0.05039 1 

181.31“ 

0.07554 1 

178.03“ 

0.07883 

i 

177.63" 

% 

”3 

_ 1 
1 


C. 02535 1 
1 

181,07“ 

0.02054 

I 

I 

178.99“ 

q /b 

1 




0.00047 

I 

159.19“ 

”4 

1 


1 



i 


% 

0.14780 ' 
1 

160.48“ 

0.14548 1 

160.42“ 

0.14364 

I 

160.30” 


0.04973 j 

181.88“ 

0.05680 1 

181.47“ 

0.05937 

I 

181.34“ 


1 

1 


0.000618 1 
1 

122.60“ 

0.000565 

I 

I 

242.35“ 

q , /b 



1 


0.000081 


- 6.87“ 


1 


1 



I 



1.0 I 

0.0* 

1.0 1 
1 

0.0“ 

1.0 

I 

0.0“ 


0.20484 j 

175.75“ 

0.18586 1 
1 

175.28“ 

0.18223 

i 

175.24“ 


1 


0.05215 j 

176.65“ 

0.05252 

I 

I 

176.50" 


~ 


- 1 

1 


0.01955 

I 

J. 

177.66“ 


TABLE 5.2 CONCLUDED 
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FIGURE 5-1 Flow Chart for Locating Eeros of (5-8) 




















Chapter VI 


AEROELASTIC MODES USING AIRLOADS FROM INCOMPRESSIBLE STRIP THEORY 

A. The Effect of Steady Drag on Flutter of a NonllftlnK Wing 

Before considering steady-state deflections due to lift, a thorough 
understanding of the stability behavior of the cantilever wing at zero 
steady lift is needed. With and at zero, the fore-and-aft 

bending degree of freedom is dynamically uncoupled from vertical bending 
and torsion motions, and the system analyzed in Ref. 1 results. In the 
zero lift case, then, solutions for stability involve 2n aeroelastic 
modes consisting of coupled motions In and ; the remaining n 

modes represent uncoupled free vibration in each of the assumed modes 
in v^ . As demonstrated In previous chapters, the flutter conditions 
obtained by this assumed mode analysis compare favorably with Ref. 1 
results over all practical combinations of the parameters M , P , , 

A , S , and C . Owing to this good agreement, the results and conclu- 
sions of Ref. 1 apply here as well, yet the assumed mode solution method 
still is useful in providing additional insight into the flutter behavior 
of the nonlifting wing. 

The parameters M , , A , and S offer no suprising effect, 

and most importantly here, the nature of their influence is not altered 
by the inclusion of drag. The Ref. 1 results Indicate that an increase 
in the elastic axis - A.C. offset given by A is destabllzlng, an 
increase in the sectional C.G - elastic axis offset given by S is 
destabilizing, an increase in the radius of gyration parameter i^ is 
stabilizing, and that the flutter speed is approximately proportional 
to the square root of the mass ratio parameter M . Since Ref. 1 

! 
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establishes that the Influence of drag is not sensitive to these four 
parameters, more detailed study in this area Is not considered here^ 

The most curious finding of Ref. 1 involves the effect of C in 
conjunction with the aspect-ratio parameter P , which is the product 
of the vertical bendlng-to-torsion stiffness ratio with the inverse 
square of the geometric aspect ratio. The effect of steady drag on 
flutter speed is stabilizing for smaller aspect ratios (larger P) but 
destabilizing for larger aspect ratio wings. The reversal of the effect 
of drag on flutter occurs near P •• 0.01 , and the behavior In this 
neighborhood, including flutter mode shapes, appears to be quite 
interesting- No conclusions regarding physical causes of this phenomenon 
were made in Ref. 1 , however. 

Figure 6-1 Is a reproduction from Ref. 1, showing the effect of 
drag on flutter speed as a function of P for intermediate values of 
M , 1^ , A , and S . Numbers in parentheses on the abscissa give the 
true aspect ratio for a typical value of the ratio EI^/GI^ =1.6 . 

Clearly wings of practical interest include the region within which the 
effect of drag on flutter appears to be most interesting. 

To help gain a better physical understanding of the behavior near 
P = 0.01 , flutter solutions for this same example have been fovind via 
the simple harmonic method of Chapter IV, over the range 0.002 ^ P ^ 0.02. 
Results appear in Fig. 6-2 that show flutter speeds and flutter mode 
shape amplitudes and phase relations as functions of P . Since a 
finite chordwise to vertical bending stiffness ratio T must be 
specified, and the effect of drag on flutter depends on T as in 
Fig. 4-15 the value x = 50 was used to allow adequately for the 
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T -»■ » behavior inherent in the Ref. 1 formulation. Three assumed 
inodes In each degree of freedom are used. 

The flutter speeds In Fig. 6-2a closely match the Fig. 6-1 results. 
The flutter mode shapes include participation by the three generalized 
displacements q , q , and ci , with the remaining assumed modes 

<Pl Wj^ -Wj 

contributing negligibly to the motions. In Fig. 6-2b the amplitude and 
phase of the two assumed bending inodes at flutter are showti for unit 
magnitude and zero phase angle of the first assumed torsion mode. 

It can be seen that for any C the sharp drop in flutter speed 
that comes with decreasing P is accompanied by a sudden change in the 
flutter mode shape. The amplitudes of the two bending modes merge, and 
the second bending mode undergoes a large phase shift. Further decrease 
in P gives a gradual separation of the assumed bending mode amplitudes 
with the fundamental mode again becoming dominant. 

To help visualize the physical appearance of these flutter mode 
shapes, phasor diagrams of the spanwlse distribution of bending displace- 
ments, for C = 0.02 , are given in Fig. 6-3 for five values of P . 
Arrows depict the modal generalized displacements from Fig. 6-2, and 
the curves give the relative displacements along the span and their 

phase referenced to q, . For both P = 0.02 and P = 0.002 all 

•^1 

stations are nearly in phase and the mode shape is dominated by the 
first assumed bending mode. But for the intermediate values of P , 
where the transition in phase of the second bending mode takes place, 
the displacements at different locations along the span can be over 90® 
out of phase. 
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The behavior of flutter nodes in this range of P offers the 
greatest discrepancy found between the results reported In Ref. 1 and 
Che assumed mode solution, suggested by Fig. 2-1. For P < 0.04 in 
Ref. 1, collocation at only five spanwise stations was used, and the 
mode shapes in bending and torsion were permitted to have spanwise phase 
differences. It was found, however, that these phase differences never 
exceeded a few degrees, in contrast to the results presented here. 
Possibly the use of only five spanwise collocation points did not allow 
enoLigh freedom to represent the flutter mode shape transition found 
using assumed modes. In any case good agreement between flutter speeds 
and frequencies is still observed for the two methods. 

Figure 6-2 seems to indicate that the second assumed bending mode 
plays a significant role in the reversal of the effect of drag on flutter 
near p = .01 , which coincides with the natural frequency of this mode 
crossing the flutter frequency. Interaction of actual aeroelastie modes 
is masked by the limitations of the solution method, however, which only 
gives neutr.ally stable solutions. In order to better understand these 
results, the Laplace transform approach of Chapter V is used to allow 
tracing ail of the aeroelastie modes in the complex plane for speeds 
from zero into the supercritical range. 

Figures 6-4 (a) - 6-4 (h) give root locus diagrams for increasing 
speed at eight representative values of aspect-ratio parameter P . 

Each locus originates for U = 0 at one of the normal modes of free 
vibration of the structure, which are easily calculated in terms of the 
uncoupled assumed bending and torsion modes. Zero drag branches are 
shown in all of the figures, with loci corresponding to C 0 added 
where their behavior differs significantly from that for zero drag. 
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The first two figures represent stubby low-aspect-ratio wings, for 
which the strip-theory aerodynamic assumption is certainly inadequate. 

Due to the manner In which the Laplace variable p is nondimensionalized, 
the predominantly first torsion normal mode of free vibration remains 
essentially fixed on the itu axis near 1.6 on all of these diagrams. 

The predominantly bending normal modes move down the iu) axis as P 
decreases, since their natural frequencies decrease relative to the 
torsion frequencies. In Fig. 6-^ (a) the bending branch of the aero- 
elastic modes leads to flutter, whereas in Fig. 6-4(b) the torsion 
branch eventually becomes unstable. The normal mode having the third 
lowest natural frequency, predominantly the second bending assumed mode, 
occurs well up the _lu axis and off these two diagrams and has negligible 
influence on flutter. These low-aspect-ratio cases show entirely two- 
degree-of- freedom behavior and closely resemble the root locus given 
fey Edwards (Ref. 11) for a typical section in plunge-and-pitch motion 
in incompressible flow. 

Figures 6-4 (c) and 6-4(d) represent values of P just above the 
condition where the effect of drag on flutter reverses. Although the 
flutter phenemenon is still similar to that for larger P , the third 
normal mode frequency has now decreased sufficiently to appear on the 
diagram, and it produces .a branch which does not lead to flutter for all 
values of drag. 

According to the flutter curve of Fig. 6-2, for P = 0.01 the 
cases C = 0.02 and C = 0.04 result in decrease of flutter speed 
from the zero-drag condition whereas for C = 0.01 it is still increased. 
Figure 6-4(e) gives the P = 0.01 root locus, which reveals that the 
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aeroelastic mode emanating from the predominantly second bending normal 
mode now plays an important role. The next illustration. Fig. 6-A(f) 
with F = 0.007, gives this branch as becoming unstable for all C . 

In the final two of these illustrations, with P = 0.005 and 0.002 , 
the flutter phenomenon appears to be returning to the type of behavior 
seen for small-aspect-ratio wings, with the second bending contribution 
assuming a lesser influence. In Fig. 6-4 (h) a fourth normal mode, the 
third bending mode, has made its appearance but does not noticeably 
influence flutter. 

The nature of the aeroelastic modes on the various branches of the 
root locus diagrams can be clarified by looking at their mode shapes. 

In fig. 6-5, phasor diagrams are used to show C = 0 mode shapes for 
each branch of the .005 locus {Fig. 6-4(g)). At selected speeds both 
subcritical and supercritical, the generalized displacements are shown 
with the phase angle of taken as zero. 

The dependence of drag's effect on flutter upon P is also given 
for a second configuration having a smaller mass ratio M = 9.4 , more 
representative of light aircraft and sailplanes. Shown in Fig. 6-6, 

Che basic behavior resembles the first configuration, with certain 
differences. For example, the reduction in flutter speed by drag for 
small P is moderate relative to that for the larger mass ratio; this 
is also seen in the tabulated Ref. 1 results. Also the transition as 
P reduces through the 0.01 region is much less severe. Again the 
reversal coincides with the crossing of the second assumed mode frequency 
below the flutter frequency. 

In conclusion, assumed mode results attribute the reversal of the 

effect of drag on flutter speed to the Interaction of the structure's 

t 

lie 


I nil I mil 







«?:-conrl naturil bending mode with the basic cantilever wing flutter 
mechanism involving the first bending and first torsion modes. The 
exact reason that drag increases the flutter speed for V above this 
apparent "resonance" condition and decreases It for smallt r P 1s not 
evident in these results, and like many .seroelastic phenomena may not 
have a simple physical explanation. It does appe.ar, liowever, that the 
drag force enhances Che coupling of the second bending mode into the 
flutter behavior and thereby magnifies Its .alre.sdy present effect on 
flutter speed. 

Since the Ref. i tabulated results give a fairly complete picture 
of flutter of the uonliCtlng case in the presence of steady drag and 
are not disputed by current results, further work herein is directed 
towards the more general case Involving steady deformations due to lift. 

b’ Effect of Steady Deformations Due to Lift 

The good agreement In prediction of dynamic stability between the 
current analysis and the Ref. 1 collocation method for nonlifting wings 
witl> a steady dv.ag force included furnishes confidence Ch.at the modal 
scheme will be successful Cor steady lifting conditions. The effect of 
steady deformations due to lift is Incorporated into the dynamic 
stability analysis by the same means as the steady drag effect-“ namely 
through coefficients of the stiffness matrix detertii.'ned in a separate 
nonlinear solution for the steady-state deflections. Thus, the .agreement 
indicates that the scheme allowing sm.all time dependent perturbations 
about a static deflection is working properly. W\en steady lifting 
deflections ati* introduced through .a nonzero root angle of attack ct , 
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Che fore-and-aft bending is no longer dynamically uncoupled. 

Coupling in both elastic terms (the stiffness matrix coefficients 
containing and 4>^) and aerodynamic coupling terms (arising from 

unsteady aerodynamic force components in the x-dlrectlon) now appears 
in the fore-and-aft banding dynamic equations. The 3n aeroelastic 
modes will consist of coupled motions In '^ 2 . » 'f'j > ^^2 ' 

Three basic vfing configurations are selected to illustrate the 
effects Introduced by steady-state lifting deformations. Parameters 
M , , A , and S are taken the same as in Fig. 6-1, and aspect ratio 

parameter P is assigned three different values in order to consider 
wings of large, moderate, and low aspect ratio. For large aspect ratio, 
P = 0.005 is chosen to provide a case for which steady drag decreases 
the flutter speed (Fig. 6-1) . A moderav-e aspect ratio example with 
P = 0.02 having an increase in flutter spet.* due to steady drag, and a 
low-aspsct-ratlo case, P = 0,1 , are also included. For the typical 
bending-to-torsion stiffness ratio 



these examples correspond to actual retangular planforms having aspect 
ratios of 17.89 , 8.94 , and 4 respectively. The bending stiffness 
ratio T now becomes important, and is given the nominal value 25. 

The essential features of the flutter behavior encountered when steady 
deformations enter can be Illustrated by using these basic configurations 
as examples. 

The flutter stability boundary for the raoderate-aspect-ratio example 
is shown in Fig. 6-7 for C = 0 and C = 0.01; the steady bending 
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displacement vr^ of the wlngtip Is the measure of steady lift. An 

alternative would be to show flutter speeds as a function of a , but 

this is a poor means for comparing curves having different steady drag 

and gives no information about the elastic steady deformations. A 

better way to indicate the steady flight condition would be the total 

lift force on the deformed wing nondlmenslonallzed, for example as 
LiJl ^ 

7 _.- ■ . As can be seen in Fig. 6-8 this dimensionless total lift para- 
4E1 0 

X 

meter, which depends on the steady twist » varies for constant 

wingtlp deflection as C changes for points along the stability 
boundary of Fig. 6-7. This is because the drag force alters the rela- 
tive tb and w distributions for the same total steady lift. But 
a o 

since this effect is small, and gives the best indication of the 

magnitude of the steady equilibrium deflections, this deflection is used 
to Indicate the steady lift condition. In Fig. 6-7, the semispan of the 
wing is about 9 semichords, and steady deformations well exceeding the 
limits of the moderate displacement beam theory are therefore shown. 

This demonstrates that flutter solutions can be found for arbitrarily 
large steady displacements and that It is a matter of practical engineer- 
ing judgment to recognize when the assumptions made in the derivation 
of the equations have been violated. 

For C “ 0 Fig. 6-7 indicates that the flutter speed reduces 
continuously with increasing steady lift until a maximum reduction of 
about 13% is achieved at an excessively large deflection of 4 semichords. 
The frequency at flutter reduces monotonlcally with w^ ; this effect 
is generally observed for all wing configurations. 



With steady drag Included, the situation for the deformed wing is 
not as simple as it is for the a " 0 case. Ae exhibited by the non- 
linear steady deformations shown in Chapter IV, drag can greatly alter 
the deformation state associated with a given speed and root angle of 
attack. Even more importantly, drag aignificantly reduces divergence 
speeds, possibly to leas than the flutter speed. This is revealed in 
the stability analysis when the nonlinear steady solution blows up 
before dynamic instabilities appear. 

In Fig, 6-7 divergence speeds found by the linear VT dotamlnant 
(2-29) are also indicated for several values of C . Results of a 
dynamic stability analysis with C * 0.01 also appear. The flutter 
speed for a = 0 and C - 0.01 is the single point on the ordinate 
and is greater than the divergence apeed for the same drag. During the 
search for a neutrally stable oscillating condition for the very small 
angle of attack a « 0.001 rad. and C • 0.01 , the nonlinear steady 
displacement solution was sensitive to U in the neighborhood of flutter 
owing to the proximity of divergence. For specified larger angles of 
attack, the C » 0.01 flutter boundary is found without difficulty. 

The slight increase in flutter speed due to drag observed for a = 0 
appears to be preserved in the presence of steady deformations. For 
decreasing Ot , however, a point is reached as a ->■ 0 at which the 
steady displacements are still nonzero at flutter. For smaller steady 
deflections divergence becomes the instability encountered for increasing 
speed. For C = 0.02, divergence was observed for all lifting conditions 
with no flutter points found. 
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As steady deformation increases from zero, the flutter mode shape 
undergoes a smooth transition, originating with the same zero lift mode 
shapes shown in Fig. 6-2 for P = .02 . For zero drag, the amplitudes 
of modal generalized displacements which contribute significantly to 
the flutter mode are shown in Fig. 6-9, normalized to the first torsion 
mode amplitude. Their phase angles relative to zero phase for q. 


^1 

appear in Fig. 6-10. Participation by the first chordwlse bending mode 
Increases steadily with Increasing , reflecting the increased strength 
of the elastic bending-torsion coupling. Vertical bending motions are 
increasingly dominated by the first assumed bending mode, and the 
contribution of the second torsion assumed mode increases significantly. 

The flutter mode shapes at one steady lifting condition are presented 
in Fig. 6-11 in a form giving a clearer physical description of the 
motion. Above the phasor diagram (which contains the information given 
in Fig. 6-9 and 6-10) is a sketch of the cyclic path traced In the 
x-z plane by wing sections at the wlngtip and at mldsemispan. Points 
where the first torsion assumed mode is at phase angles of 0° , 90“ , 

180“ , and 270“ are located. These diagrams emphasize the three-degree- 
of-freedom nature which flutter can have when steady deflections are 
present. 

The low-aspect-ratio example (P = 0.1 and /R. ~ 4) shows only minor 
effects upon its flutter characteristics due to steady deformations, as 
might be anticipated. As given in Fig. 6-12, even for the extreme 
condition a = .12 rad. yielding a 1.7 semichord tip deflection at 
flutter, there is only a 0.62% reduction in flutter speed due to lift. 
The flutter mode shapes (Fig. 6-12) undergo little change, with a slight 
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contribution by q the only new feature. Owing to these unremarkable 
results, further work centers on the inoderate-and-high-aspect-ratlo 
examples, where the effects of steady deformationr: and chordwise forces 
are significant and Interesting. 

For the high-aspect-ratio example with P = 0.005 , the flutter 
behavior Is quite different when steady deformations enter. Figure 
6-13 shows the dependence of flutter speed upon the steady wingtip 
bending deflection for two drag cases, C “ 0 and 0.01 . The minimum 
flutter speed in this case is over 20% below its undeflected counterpart, 
but quite Interestingly as w^ 0 the stability boundaries do not 
converge continuously to their respective zero lift flutter speeds but 
approach lower points on the ordinate. 

The flutter mode shapes. Figs. 6-14, 6-15, explain this new 
behavior. As the steady deflection becomes small, the flutter mode 
becomes dominated by the first chordwise bending mode, and as « 0 

this type of Instability approaches simple free uncoupled vibration in 
this degree of freedom. For moderately large steady deflection the 
flutter mode shape closely resembles that for the medium-aspect-ratio 
wing. Figs. 6-9 and 6-10. 


Flutter mode shapes for t\io steady-lift conditions are diagrammed 
in Fig. 6-16 using the same technique as in Fig. 6-11. Relative to the 
P = 0.02 example, this wing shows a greater amount of participation in 
fore-and-aft bending, and a greater contribution from the second vertical 
bending mode. This mode's contribution causes the noticeable difference 
in the eccentricities of the elliptical paths traced in the x-z plane 


by different wing stations. The phase relationship among 
and q appears similar In this and other examples. 
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The most noticeable effect of steady drag again is Its reduction 
of the divergence speed. Figure 6-17 gives divergence speeds found by 
the linear determinant (2-29) together with zero lift flutter speeds 
for Increasing steady drag for the moderate- and high-aspect-ratio 
cases. Divergence is clearly more important relative to flutter for 
the larger aspect ratio. 

The effect of steady drag on the dynamic behavior (Fig. 6-13) 
appears to cause only slight adjustments to the zero drag flutter results 
for any steady deflection. The flutter mode shape amplitudes and phase 
angles show but small changes for the rather large drag C = 0.01 . 

The flutter speeds fall just below the C = 0.01 divergence speed. 

The stability boundaries In Fig. 6-13 do not allow for structional 
damping and give no feel for the degree of stability at speeds near 
flutter. To gain a better understanding of the type of instability 
that has been found with steady deformations, and also define the over- 
all aeroelastlc behavior, the Laplace transform approach detailed in 
Chapter V is applied to the large- and moderate-aspect-ratio examples. 

The true aeroelastlc modes are then conveniently traced in the complex 
plane at any flight conditions using root locus diagrams. 

Before showing the root-locus results, it is interesting to see 
how stability is suggested by the simple harmonic method with nonzero 
structural damping assumed. The C = 0 stability boundary of Fig. 6-13 
is reproduced in Fig. 6-18, to which stability boundaries for three 
values of structural damping are added. The sizeable increase in 
flutter speed with g hints that the instability is not of severe 
nature, and as w^ becomes small the predominantly fore-and-aft bending 
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motion is obviously very lightly damped. For small deflections with 
structural damping included, there is a change in the flutter mode back 
to the basic bending-torsion type encountered for zero steady lift. In 
spite of structural damping's stabilizing influence, there still exists 
the possibility of a reduction in flutter speed at high load factor 
flight conditions. 

Root locus diagrams depicting the dynamic stability of the moderate- 
and large-aspect-ratio-wing examples with steady deflections appear in 
Figs. 6-19 and 6-20 respectively. Dashed lines in these figures show 
the zero steady lift loci of roots for increasing speed, which are the 
same as the zero drag diagrams of Figs. 6-4 (c) and 6-4 (g). Solid lines 
trace the elastic modes for selected constant speeds as the angle of 
attack a is varied, and originate for a = 0 at nonlifting roots 
corresponding to these speeds. 

In addition to the nojnnal modes of free vibration involving vertical 

bending and torsion, previously seen on the lu axis in Figs. 6-4, the 

first fore-and-aft bending mode natural frequency now also must be 

included. This normal mode remains an uncoupled, undamped aeroelastlc 

mode at all speeds at a = 0 , but a family of constant speed branches 

emanates from this root with steady lift included. Since Its natural 

frequency is a factor /r larger than the first vertical bending 

assumed mode frequency u , this new normal mode lies on the iu) 

”l 

axis approximately at a multiple /F of the lowest normal mode frequency. 
This falls above the zero-lift flutter frequency for the moderate 
aspect ratio wing, but less than it in the high-aspect-ratlo example. 
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A general feature of Fig. 6-19 and 6-20 Is that as steady lift 
Increases from aero, one group of constant-speed branches tends to 
stabilize and rapidly increase In frequency, whereas another lower 
frequency family of root paths reduces In stability and decreases in 
frequency. This latter group is responsible for the stability boundaries 
of Figs. 6-7 and 6-13. In the high-aspect-ratio example, they originate 
from the lowest chordwise bending normal mode, but for the moderate- 
aspect-ratio wing they originate on the zero-lift torsion branch of the 
root locus. The constant-speed branches which originate on the first 
and second vertical bending zero lift paths, for both examples, do not 
show much sensitivity to steady deflections. 

To illustrate more clearly the role played by the chordwise bending 
upon stablltcy, the large-aspact-ratlo example is modified by increasing 
the bending-stiffness ratio t from 25 to 60, This raises the dimension- 
less natural frequency of the first chordwise bending mode from 0.6216 
to 0.9629, which is greater than the zero lift bending- tors ion flutter 
frequency of 0.8921. The stability boundary calculated for this modified 
example appears in Fig. 6-21 and the associated flutter modes are 
presented in Figs. 6-22 and 6-23. The root locus obtained via the 
Laplace transform method appears in Fig. 6-24 (a) and the true stability 
of constant-speed branches yielding instability is better depicted in 
6-24 (b) using the damping ratio C • 

The discontinuity in the stability boundary is only a consequence 
of the solution procedure of Chapter IV and covers a region where solu- 
tions that do exist could not be determined. It is due to an interaction 
of the predominantly second bending, stable aeroelastic mode with the 
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roots that yield instability. TVie simple harmonic solution method 
involves fixing the angle of attack o and searching for the neutrally 

I* 

stable lifting condition whose speed matches the calculated flutter 
speed. In Fig. 6-21, the gap in the curve falls between points found 
for a = .0095 and .01 rad. 

For all points on the left segment of the stability boundary, the 
simple harmonic solution yielded an additional highly damped eigenvalue 
whose frequency was slightly below the flutter frequency. On the right 
portion its frequency was above the flutter frequency. From the root- 
locus diagram (Fig. 6-24 (a)) this highly damped eigenvalue can be 
identified as the predominantly second vertical bending aeroelastic 
mode, and the discontinuous behavior In Fig. 6-21 coincides with the 
crossing of frequencies as the downward moving constant speed loci 
associated with flutter pass the stable second bending aeroelastic mode 
frequency. By correlation with the section VI-A discussion of the effect 
of P and C on flutter, the general effect of the second bending mode 
appears to be destabilizing when its frequency is just below that of 
flutter, and stabilizing when its frequency is just above. 

The flutter mode shapes (Fig. 6-22, 6"23) reflect an interesting 
transition as flutter frequency drops below the second bending aero- 
elastic mode frequency for increasing w^ . To the left of Che 
discontinuity the mode shapes resemble the zero-lift large-aspect-ratio 
flutter behavior, to which they converge as a -*■ 0 . For larger , 
though, the modes closely resemble flutter mode shapes for the t = 25 
case (Figs. 6-14, 6-15) which in turn resemble mode shapes for the 
TOoderate-aspect-ratio example. This sliularity is most evident in q 

w„ 



In Fig. 6-25 the flutter mode shapes are dlagranmed for two lifting 
conditions to assist physical visualization of the motion. The larger 
steady deflection illustration is similar to those given in Fig. 6-16. 

But the smaller steady deformation mode shape, with its large contribu- 
tion by the second assumed bending mode, has the wingtip rotating in a 
clockwise direction about its elliptical path for a counterclockwise 
motion at midseraispan. In spite of the relatively small steady deflec- 
tion (a deflection of 1.226 semichords at the wingtip for a semispan of 
roughly 19 semichords) the motion Is quite three-dimensional, indicating 
a significant inertial contribution to flutter in fore-and-aft bending. 

The three root-locus diagrams {Figs. 6-19, 6-20, and 6-2t) exhibit 
die basic effect which steady deflections have upon dynamic stability 
when incompressible strip-theory airloads are used. The basic zero- 
steady-lift bending-torsion flutter root together with the first fore- 
and-aft assumed mode produce a pair of constant-speed branches v one of 
which rapidly stabilizes and increases in frequency while the other 
decreases In frequency and becomes unstable for speeds below the a = 0 
flutter speed. These latter aeroelastlc modes are lightly damped, and 
generally the onset of flutter at constant speed for increasing steady 
deflections would not be as severe as for that encountered with increasing 
Speed. The reduction of flutter speed with steady deformations has been 
noted to be greatest when the first fore-and-aft natural frequency Is 
near to the basic bending- torsion flutter fTequency. 

A better comprehension of the various aeroelastlc modes comes from 
inspecting mode shapes at both subcritical and supercritical conditions. 

In Fig. 6-26, mode shapes for selected points along the U = 6 and 7 
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branches in FIr. 6-19, the moderate-aspect-ratio example, are shown in 

phasor form. Likewise for the high-aspect-ratio wing having T = 60 , 

Fig. fi-27 presents mode shapes for certain steady deflections along the 

U = 6 and 7 paths from Fig. 6-24. These mode shapes disclose that 

the branches which stabilize and increase in frequency consist largely 

of motion in the first assumed torsion mode, having a small q 

contribution locked in a characteristic phase relationship with q. of 

^1 

near 180“. The branches producing instabilities, furthermore, have q 

^1 

nearly in phase with q^ and give the previously shown flutter mode 

■ 1 

shapes as they cross the iio axis. Made shapes (Fig. 6-27) for the 
essentially second bending aeroelastic modes show the dominance of the 
second assumed vertical bending mode in this branch, which remain nearly 
fixed in Che complex plane at a frequency close to the second assumed 
bending mode natural frequency. 

The role of the bending-stiffness ratio T in flutter of lifting 
wings is next examined. Figures 6-28 and 6-29 give flutter speeds and 
frequencies found for the same moderate-aspect-ratio example used earlier 
compared with solutions for different values of T . For t = 12 the 
dimensionless natural frequency of the first chordwise bending mode is 
0.86124, which is only slightly higher than the zero-lift flutter 
frequency of 0.8332, and its stability boundary shows the most marked 
decrease in flutter speed as w^ increases from zero. The other 
extreme, T = 1000 , has a dimensionless natural frequency in chordwise 
bending of 7.862, yet some decrease in flutter speed still occurs. 

Elastic deformation in bending about the major principal axis of Che 
airfoil cross section should be virtually suppressed and all fore-and-aft 



motion connected with vertical bendlng'-torsion elastic coupling. Calcu- 
lations are impractical for larger T since numerical problems begin 
to appear, such as convergence difficulties with the nonlinear steady 
deflection iterative solution. 

The uniform decrease in flutter frequencies as T is decreased 
reflects the additional inertia from the larger fore-and-aft motions 
brought about by the reduced elastic stiffness in bending about the 
airfoil major principal axis. Decreasing flutter frequencies for 
increasing steady deformations likewise should be due in part to the 
increasing contribution of q , although the increase in the relative 
participation of q^^ is also a factor. 

In Fig. 6-30 the stability boundaries shown earlier for the large- 
aspect-ratio example with X = 25 and 60 are compared along with 
curves having T = 10, 40, and 200. The appearance is complicated by 
the interaction of the flutter modes with the predominantly second 
vertical bending aeroelastic mode as discussed earlier, and by fore-and- 
aft bending natural frequencies sufficiently low to cause convergence 
as a -*■ 0 to speeds below the true nonlifting flutter speed. 

For both T = 200 and 60 the flutter frequencies descend through 
Che range of the second vertical bending mode frequency as is 

increased, causing the discontinuous stability boundarie.s. The flutter 
frequencies (Fig. 6-31) reflect the role of this second bending mode, 
whose dimensionless natural frequency is 0.77904 and (as already shown) 
a stable aeroelastic mode with approximately this frequency exists at 
flutter. As mentioned in discussing the T = 60 results, these curves 
are actually continuous, but the simple harmonic solution procedure 
could not produce results within the gaps. 
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For the t = 40 , 25 , and 10 cases, stability boundaries converge 

to free vibration in the first fore-and-aft mode as a -* 0 ; their 

frequency curves converge to the respective natural frequencies of this 

mode. In the case T = 40 , this frequency is 0.7862, quite close to 

the second:, vertic,al bending assumed mode natural frequency 0. 77904. 

This near resonance caused difficulty in precisely locating the neutral 

stability curves for small , ns the damping became exCremelv small: 

of course this trouble is uniraport,int since even small structural damping 

would raise flutter speeds considerably here. 

The stiffness in chordwisc bending is thus a factor in the flutter 

of lifting cantilever wings. Although stability is most adversely 

affected for uj near the zero-lift flutter frequency, this effect 
'’l 

appears over a wide range of T . Hlgher-,aspecC-raclo-vfings apparently 

experience a greater decrease in flutter speed with steady deformation.';, 

given that te is sufficiently near the n = 0 flutter frequency. 

'"l 

The mass ratio M is the only parameter which can change for a 
specific wing, as it depends upc n altitude. In Fig. 6-32 the dependence 
of both flutter and divergence speed.s on M for the moderate-aspect- 
ratio example appears for the complete range of mass ratios of practical 
interest. Divergence speeds for several values of steady drag together 
with zero-lift flutter speeds for C = 0 and 0.01 are compared. In 
addition zero drag flutter speeds for steady lift giving 


V 
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are shown for M = 10 , 40 , and 100 . This same Infortnat ion Is 
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depicted in Fig. 6-33 for the large-aspect-ratio example except wlngtip 
deflections of 2 and 6 seraichords are used. 

Divergence speeds for all drag values vary exactly as , as can 
be seen, from the divergence determinant (2-29), where M appears only 
as a product with the inverse square of divergence speed. The C = 0 
divergence speeds for these two examples are Identical, but decreases 

more rapidly with C for the higher aspect ratio. 

For mass ratios above about .5, nonlifting flutter speeds behave 
.approximately as , A minimum in flutter speed Is found near M = 5 . 

whith a rapid asymptotic rise to infinity following a further decrease 
in M . Thi.s reflects well known results for incompressible flow (Ref. 5, 
page 247). A practical application for mass ratios sufficiently small 
to be theoretically free from flutter is the stability of hydrufoil.s 
used in high-speed marine transportation. Divergence would be the type 
of instability encountered. This conclusion emph.asizes the lrapovt.ancc 
of allowing for steady drag since its effect on divergence is propor- 
tionally the same for any mas.s r.itlo nnd the decrease of divergence speed 
with C can be considerable. 

Divergence Is more important relficive to flutter for mass ratios 
around 10 than for higher mass r.atlos; this Indlc.ates a greater likelihood 
of light aircraft .and sailplanes experiencing divergence. Calcnlntion 
of divergence speeds for a high-performance sailplane appears ns n likely 
application In which allowance for steady drag effects would be deslreahle. 

The effect of steady deformation on flutter .appears insonsirive to 
M , For equal tip deflections, reduce.s by about the same proportion 

at each of the three mass ratios in the two figure.s. For the high-.aspoct- 
ratio wing, these deflected flutter speeds correspond to the neroelastlr 
mode connected with chordwlse bending, ns discussed for M - 40 . 
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The effect of changes In Che remaining parameters A , S , and 
has also been investigated for the steady deformed case, but no new 
interesting behavior was found. Their influence at zero lift, as 
tabulated in Ref. 1, Is carried over into the deflected flutter behavior. 
Stability boundaries and flutter mode shapes for increasing steady lift 
show Che same alterations of the basic flutter phenomenon as is demon- 
strated by Che examples used here. 

To clarify Che strong influence which steady chordwise loads can 
have upon divergence, particularly for high aspect ratios, the dependence 
of upon C as aspect ratio parameter P is varied appears in 

Fig. 6-34. The corresponding zero-lift flutter speeds for this 
configuration (Fig. 6-1) can be compared directly. Unquestionably, for 
higher aspect ratios (P < 0.01) drag forces typically encountered in 
flight (i.e., C = 0.0025) can cause divergence to become no less 
important for flight safety than primary bending-torsion flutter. 

C. Effect of Unsteady Leading-Edge Suction Forces from Two Dimensional 
Incompressible Flow 

Substituting the terms in (5-20) for their counterparts in the aero- 
dynamic matrix (5-3) allows inclusion of the linearized unsteady chord- 
wise forces arising from incompressible potential strip theory as 
described In Chapter V. This effect is present only for oi 0 and 
will grow as angle of attack is increased. The initial trend of stability 
as a increases from zero should be unchanged from the results already 
discussed, but at higher steady lifting conditions leading edge suction 
should have a noticeable influence. 
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Figure 6-35 displays how allowance for the linearized unsteady 
propulsive force affects the destabilizing, constant speed. Increasing 
steady lift branches of the root locus In Fig. 6-19, which 1s the 
moderate-aspect-ratio-wlng example. For this case a tip deflection of 
just one semlchord is a very large deformation assuming a conventional 
bending-to-torsion stiffness ratio EI^/GI^ , since the semispan would 
be about 9 seraichords. Thus larger deflections shown are only of 
academic interest, since they exceed in practice the limits imposed by 
several underlying assumptions. 

Figure 6-36 recasts the Fig. 6-35 information into a better form 
for inferring stability; the damping ratio of the same root branches is 
plotted for wingtip deflection. Statilllcy is slightly increased when 
leading-edge suction is added, but for practical deformations its effect 
is really not too significant. A stability characteristic of deformed 
wings, also evident In the format of Fig. 6-36, is Che reduced severity 
of flutter onset with increasing speed at constant for higher load 
factors. 

Figure 6-37 gives the manner in which the root locus for the high- 
aspect-ratio example (Figure 6-20) is altered by allowing for linearized 
unsteady suction forces. Shown is the destabilizing family of constant- 
speed branches associated with the first fore-and-aft bending normal 
mode. Again true stability with steady wingtip deflection Is also 
shown (Fig. 6-38). Also added are the constant-speed branches which 
stabilize and Increase in frequency as steady deformations Increase. 

Once again the unsteady suction force has a slightly stabilizing 
influence. Given the same ratio EI^/GI^ , this example has twice the 
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span in semlthorda of that for the moderate-aspect-ratio example, hence 
las'ger wingtlp deflections in semichords can be tolerated. After 
allowing for this factor, the unsteady chordwlse potential forces do 
not appear to have any greater Influence on stability for the larger 
aspect ratio. 

Figure 6-39 shows how stability of the modified hlgh-aspect-ratio 
example (having T = 60 , presented In Fig. 6-24) depends on steady 
deformations and on the linearized unsteady propulsive force. Comparison 
with Fig. 6-38 shows a similar effect due to leading-edge suction, but 
steady deformations definitely have a stronger and more Immediate 
destabilizing effect for the stlffer chordwise bending case. This is 
due to the proximity of the fore-and-aft bending normal mode frequency 
to the zero-lift bending-torsion flutter frequency. 

Inclusion of the linearized unsteady leading-edge suction terms 
derived In Section T) of Chapter V is thus found to be stabilizing for 
deformed wings. Aeroelastic modes involving a primary contribution by 
the first chordwise bending mode appear to be most affected by these 
terms. With their inclusion, potential flow strip theory has been 
fulJ.y exploited for this problem, and further improvement in the aero- 
dynamics involves the compressible three-dimensional loads of Chapter 
VII. 


D. Two Practical Examples 

The manner in which chordwise forces and steady deformations 
influence the aeroelastic stability of cantilever wings has thus far 
been illustrated by means of idealized examples. After one has identi- 
fied the fundamental effects, it is instructive to apply the same 
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solution techniques to typical designs. Sailplanes, having large 
aspect ratios, low mass ratio, and low operating speeds (which perirlt 
the use of Incompressible aerodynamics), are a logical choice for 
practical examples. Accordingly, two existing sailplanes are modeled 
within the approximations Imposed by the assumption of uniform spanwise 
mass and stiffness properties and strip theory airloads. Their stability 
is studied using the same techniques applied before. 

The first example Is modeled from information given about the 
Slingsby Dart 17R In Ref. 15, as summarized in Table 6.1, The approach 
taken In determining the wing stiffnesses, in lieu of looking at wing 
construction details, was to use photographs in Ref. 15, from which tip 
deflection and twist at a load factor of 4 could be measured. The 
stiffness EI^ so obtained is clearly larger than would be expected 
for the weight and type of construction; this inaccuracy is due to the 
uniform stiffness restriction, the rectangular planform, and strip theory 
steady airloads. The torsion stiffness GI^ is more reasonable. The 
ratio 



turns out impractlcally large. This number leads to a larger value of 
P than would be expected, which causes bending natural frequencies to 
be larger relative to torsion natural frequencies than would probably 
occur on the actual vehicle. It can be added that, had GI^ been 
Increased to compensate for the seemingly excessive , then the 

flutter speeds would have been absurdly high. Since wing construction 
details were not available, typical, values for the parameters A S, 
1^ , and T Were simply assigned. 
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A different philosophy was applied when modelling the second 
vehicle, the Gemini Cwo-place high performance sailplane detailed in 
Ref. 16. With construction details of the metal wing structure described 
in Ref. 16, it was possible to assign reasonably accurate values to all 
six of the dimensionless parameters from a diagram of the average wing 
cross section (Fig. 6-40). The pertinent details of modeling the 
Gemini are listed In Table 6-2. 

One notable result is a much lower bending stiffness ET^ than 
for the Dart 17R. This result will compromise the static deflection- 
load factor relationship but should favor the dynamic modeling. The 
considerable difference in bending stiffnesses between the two examples 
is revealed In Fig. 6-41, which shows how the true load factor (found 
using the respective aircraft gross weights) varies with steady vertical 
tip deflection at speeds near the expected flutter speeds. The poor 
static modeling of the Gemini should not adversely affect Its flutter 
results, however, tip deflection rather than load factor should be used 
to measure the amount of steady lift. 

The mass per unit span of the Gemini wing was taken to be the 
average for the outer two-thirds of the semispan, so as better to model 
it dynamically. The mass ratio — over twice that for the Dart 17R — 
reflects the heavier construction needed for a two-place vehicle and 
the smaller average semichord. This design (Ref. 16) intentionally 
has a higher wing loading to optimze the glide slope, with thermalling 
performance improved using full-span flaps. The chordwise bending 
stiffness was conservatively estimated yet still yielded a quite low 
value of T = 5 . The low EI^/GI^ ratio, together with Che very large 
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aspect ratio, result In the extremely small P = 0.00128 . This wing 
thus is an extreme case in the context of the examples studied earlier. 

The effect of steady drag on divergence speeds and nonlifting 
flutter speeds is shown in Fig. 6-42 for the Dart 17R example and In 
Fig. 6-43 for the Gemini. For both, the divergence speed drops below 
the zero lift flutter speed at drag values which can be realistically 
expected in flight. The excessively large P of the Dart 17R model 
moderates the drop in its divergence speed due to drag. For flutter 
only a slight dependence of speed on C is found. On the Gemini, in 
fact, the C = 0.02 flutter speed is only 0.28% less that for C = 0 . 
Its larger S tends to lower the flutter speed relative to divergence. 

Figure 6-44 depicts the stability boundary as affected by steady 
wingtlp deflection for the Dart 17R. The associated flutter mode shape 
aiiiplltudes and phase relations appear in Figs. 6-45(a), (b). The 
coupling effect between the first fore-and-aft bending mode and the 
primary zero-lift bending-torsion flutter mechanism appears to be 
responsible for a decrease in flutter speed with lift, as in earlier 
examples. The corresponding root locus diagram (Fig. 6-46) reveals that 
the first chordwise bending natural frequency by coincidence happened 
to fall almost exactly on the zero lift flutter frequency. This 
phenomenon causes the pronounced drop in flutter speed for small deflec- 
tions seen in Fig. 6-44. Inclusion of the unsteady propulsive force is 
slightly stabilizing at large deflections, as in previous cases. 

A parallel analysis of the Gemini examp3.e, given by Figs. 6-47, 
6-48, and 6-49, yields a very different response to steady deflections. 
Flutter speed decreases only slightly at representative tip displace- 
ments, and the mode shape contributions from 9,. • <1. > ^nd q show 

12 ‘’’l 
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almost no sensitivity to . The chordwlse bending contribution 
suggests the cause, made clear In the root locus diagram. The combina- 
tion of a very small P and small T produces a fundamental fore-and-aft 
bending natural frequency which is much smaller relative to the zero 
lift flutter frequency chan in any example yet treated. Furthermore, 

Che second chordwise frequency drops to near the flutter frequency. 

This second mode (rather than the fundamental) participates in flutter 
of the deformed wing, yet it Is not as strongly coupled elastically with 
the first vertical bending and torsion modes. The aeroelas tic-mode 
branches associated with the first chordwise bending mods no longer 
play an important role in stability with steady deformations. Indeed 
it has degenerated into a virtually uncoupled, neutrally damped fore-and- 
aft vibration. The linearized effect of unsteady leading-edge suction 
on this type of flutter is negligible. 

Evidently, for extreme cases such as the Gemini model with a fore- 
and-aft bending frequency much lower than the zero-lift bending-torsion 
flutter frequency, the lowest chordwise mode will not participate strongly 
in flutter. The low-frequency root branches associated with it give no 
cause for concern. Before this calculation Is accepted as a definitive 
flutter analysis of an existing sailplane, however, the cantilever root 
boundary conditions must especially be recalled. While these results 
do offer insight into symmetric flutter of the actual vehicle, the 
possibility of anti-symmetric motions involving fuselage roll and yaw 
is entirely suppressed. For the Gemini in particular, there remains 
the likelihood that anti-symmetric vertical and fore-and-aft banding 
modes may couple to produce a different type of flutter — possibly one 
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with a greater sensitivity to steady deformations. Although the results 
here are Insensitive to both steady deflections and to chordwlse forces. 
Steady drag unquesfluaably plays a critical part in divergence. 
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DART 17R SPECIFICATIOlfS (REF. 15): 


SPAN 

AREA 

ASPECT RATIO 
GROSS WEIGHT 


55.8 ft. 
149 

20.9 
800 lb. 


ASSUMPTIONS : 


WING WEIGHT 
AVG. SEMICHORD 


200 lb => m = 0.1242 

L U • 

1.33 ft. 


FROM REF, 15 PHOTOGRAPHS: 



at 4g 

=> El 


at 4g 

»> G.l 

( 


FROM THIS INFORMATION ONE CAN SPECIFY 

P^O-Ol M~9.4 (sea level) 

ASSUME THE REMAINING PARAMETERS 

= 0.25 S » 0.1 

A - 0.1 T = 25. 

EXPRESS SPEED V IN FT. /SEC. IN TERMS OF U : 

£ / T 

VT _ W tT _ rt rr 


U = ^ 

b /gI, 


V = 0.008432 V 


LOAD FACTOR AT SEA LEVEL: 


4nDV=^b«. ? 2 1 

N = (O 4- E ■= - z ~ ■ - q° ) 

W n (2j-l) 
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GEMINI SPECIFICATIONS (REF. 16) 


ASPECT RATIO 

GROSS WEIGHT 

TOTAL WING WEIGHT 

WEIGHT OF OUTER 
20 FT. 


60.5 ft. 
124 


1065 lb. 
400 lb. 


110 lb. => m = 0.1708 


ASSUME b = 1.00 ft. 

PROPERTIES ESTIMATED FROM THE TYPICAL SECTION (Fig. 6-40) 
El 444,900 Ib-ft^ 

Gl^ 402,400 Ib-ft^ 

EI^ 2,179,000 Ib-ft^ 

J 0.04697 slug-ft. 


0,2322 ft. 


RESULTING DIMENSIONLESS P.ARAMETERS 


M = 22.9 
P = 0-00128 

ia“ 


A = 0.1 
S = 0.23 
T = 5 


EXPRESS SPEED V IN FT. /SEC. IN TERMS OF U 


U = F /S- = 0*0^0335 V 

U 
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FIGURK 6-2 (a) Flutter Speeds, Frequencies, and Mode Shapes for 

Different C as P Varies, Modal Solution of the 
Fig. 6-1 Configuration 


145 




AMPLITUD 
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FIGURE 6-3 Phaaor Representation of Flutter Mode Bending Displacements 
From Fig. 6-2 for C = 0.02 at Selected Values of P 
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-1.0 -.8 -.6 -.4 -.2 0. .2 .4 

Re(p) 

FIGURE 6-4 (a) Locus of Roots as Speed Increases From Zero to Supercritical 
for Eight Values of P , With M = 40, 1 = 0.25, A = 0.1, 

S = 0.1, n = 3 
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Figure 6-4 (b) 
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P= .015 

C Up 

D 0, 7.5001 

.Oi 7.7791 

□ .02 8.0866 

0 .04 8.7146 
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FIGURE 6-4 (g) 
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Re(p) 


FIGURE 6-4 <h) 
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SECOND 

bending 

BRANCH 




FIGURE 6-5 Mode Shapes for Selected Roots on Each Branch of Fig. 6-4 (g) 
for C " 0 
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FIGURE 6-6(a) Flutter Speeds, Frequencies, and Mode Shapes for Different 
C as P Varies, Fig. 6-2 Configuration with Smaller 
Mass Ratio 
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FIGURE 6-6 (b) 
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FIGURE 6-9 Flutter Mode Shape Amplitudes for Unit |q. ■ , for the 
Fig. 6-7 Stability Boundary With C = 0 
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FIGURE 6-10 Flutter 
for the 
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FIGURE 6-11 Physical Appearance of Flutter Mode of Deflected Wing 
Fig. 6-7 Configuration With C = 0 
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FIGURE 6-12 Flutter Speed and Mode Shape Amplitudes, for Unit llx | i 
As Affected By Steady Deflections; Low-Aspect-Ratio "^1 
Example 
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ter Speeds as Affected By Steady Deformations; Large-Aspec 










FIGURE 6-17 Comparison of Drag Effect Upon Divergence and Flutter At Zero Steady Lift for the 
Moderate- and High-Aspect-Ratio Examples 
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FIGURE 6-21 Flutter Speeds As Affected By Steady Deformation, Large-Aspect-Ratio Example with 
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FIGURE 6-22 Flutter Mode Shape Anplitudes for Unit [q | Corresponding 
to Fig. 6-21 ‘^1 
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FORE-AND-AFT 
BENDING 
BRANCH 
U = 8. 



FIGURE 6-26 Mode Shapes for Selected Roots from Fig. 6-19 
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FORE-AND-AFT 
BENDING 
BRANCH 
U = 6. 


(flutter mode) 



0 . 




2.57 



q^, 


5.29 


TORSION 
BRANCH 
U — 6. 




FIRST 

VERTICAL 

BENDING 

BRANCH 

U = 6. 


FIGURE 6-27(a) Mode Shapes for Selected Roots from Fig. 6-24 
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SECOND VERTICAL 
BENDING BRANCH 
U = 7. 
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FORE-AND-AFT 
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BRANCH 
U = 7. 
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(flutter mode) 


FIGURE 6-27 (b) 
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FIGURE 6-28 Effect of T on Flutter Speeds, Moderate-Aspect-Ratio Example 



1000, 
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FIGURn 6-29 Ef 




FIGURE 6-30 Effect of t On Flutter Speeds, Large-Aspect-Ratio Example 
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FIGURE 6-31 Effect of T On Flutter Frequenclefl, Lar^e-Anpect-RAtlo Exanple 






FIGURE 6-33 Efteot of Mass Ratio on Divergentie and Flutter, targe- 
Aspett-Ratio Example 
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FIGURE 6-36 Damjilng Ratio of the Branches in Fig. 6-34 



Im (p) 


FIGURE 6-37 



Re (p) 


Effect of Allowing for Unsteady Leading Edge Suction Upon 
the Constant Speed Branches of Fig. 20 
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LOAD FACTOR 



FIGURE 6-41 Comparison of Vertical Bending Deflections of the Two 
Sailplane Wing Models for Given Steady Load Factor 
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FIGURE 6-44 Stability Boundary For the Bart 17R Model 






FIGURE 6-45 (a) 


Dart 1?R Flutter Mode Shape Amplitudes. Normalized 

“ ivi 
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FIGURE 6-45 <b) 


Dart L7R Flatter Mode 
of q. 


Shape Phase Angles for Zero Phase 
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Chapter VII 


FLUTTER VELOCITY USING AIRLOADS FROM THKEE-DlMENSIOK?AL 
SUBSONIC AERODYNAMIC THEORY 

A. Inclusion of Three-Dimensional Airloads 

The influence of steady deformations and chordwlse forces upon 
dynamic atahility of the uniform cantilever wing has been examined in 
Chapter VI using lifting airloads predicted by Incompressible steady 
and unsteady strip-theory . This approximate modeling 
of the aerodynamic loads made possible their convenient numerical compu- 
tation for any convergent, neutral, or divergent oscillations of interest. 
As a result, the iterative solution schemes of Chapters IV and V could 
be developed and a variety of wing configurations could be analyzed 
efficiently. 

The accuracy with which the incompressible strip theory results 
approximate the three-diTnenslorial compressible flow situation is next 
explored by extension to subsonic three-dimensional lifting airloads. 

Tlie simple harmonic flutter solution method described in Chapter IV is 
accordingly modified to use subsonic three-dimensional steady and 
oscillating unsteady aerodynamic loads calculated separately by the 
computer program written by Rowe, et al. (Ref. 18). Results are then 
found which demonstrate that the phenomena discussed in Chapter VI 
still occur after three-dimensional aerodynamics are Introduced. The 
role of unsteady potential chordwlse loads in flutter is Investigated 
as is ths effect of coTnprassibllity. Since the use of externally 
computed air loads requires a more cumbersome solution method, only 
enough results are sought to provide direct comparison with the Incompres- 
sible strip-theory calculations. 
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The Rowe computer program solves the pressure-downwaah Integral 
equation (Ref. 24) for compressible flow about a steady or oscillating 
planfonn. The kinematic downwash boundary condition for each structural 
mode Is enforced by collocation, at a set of user-specified points, of 
the downwash distributions associated with an assumed series of pressure 
functions. Ror this application seven collocation chords are specified 
having five collocation points per chord. Six elastic structural modes 
are input — three in vertical bending and three in torsion. The three 
torsion modes are defined for A = .1 (lateral displacements depend on 
elastic axis location) , A thorough description of the theoretical 
aspects of the subsonic kernel function program is provided In Ref. 18, 
and the programming details are documented in Ref. 19. Its capability 
for modeling trailing- and leading-edge control surfaces is not required 
in the present application. 

The unsteady potential chordwlse forces can be deduced from program 
output, as follows. For each structural mode and its downwash the 
program calculates the complex amplitude of an associated distributed 
lifting pressure difference on the rectangular planform, which in 
dimensionless form for the jth mode is, per unit q^ having steady or 
simple harmonic time dependence, 

AP.(x,y) aJfr-rr~ ^ M , . , . . , 

(7-1) ACp (x,y) = E E f (y)g^'‘\x) 


Here a^^^ are coefficient multipliers of the series expansion of 
pressure on the planform, which can be listed in the program output. 
The assumed spanvlse pressure distributions are 
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(7-2) 




sin(2v-l)<t> 

sln4 


V - 1,2,3, .. .N 


(7-3) <t> “ cos ^(^) 

N = number of downwash chords on semispan 

The assumed chordwise pressure distributions, dependent upon x only 
for a rectangular planform, are 

(7-4) g^^'fx) = / cot |- U = 1 

I sin(y-l)6 p = 2,3,4, .. .M 

(7-5) 0 = cos"^(- h 

D 

M = number of downwash points on a downwash chord 

where x is measured aft from the midchord. 

The resultant chordwise component of potential airloads, acting 
in the positive x direction, can be expressed as 

(7-6) Dp(y;t) = - ^spV* ACp(x,y;t) (x,y;t)dx - Fg(y:t) 

The first terra represents the x-component of the force, which is normal 
to the deflected chord, second terra contains Che contribution of leading- 
edge suction. It is an idealization of linearized theory, which is 
supposed to approximate the actual effects of low pressure acting around 
a curved leading edge. 

Steady and unsteady parts of the pressure and mean-surface chordwise 
slope can be separated: 
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(7-7) 


(7-8) 


ACp(x,y;t) = 


AC^ (x,y) + 


ACp (x,y;t) 


3z 

-^(x,y;t) = - 4)^(y) - 


The first term of (7-6) will then contain a steady part, a linear 
unsteady part, and a nonlinear part, as follows j 


(7-9) Dp(y;t) = HpV== dx 


+ (AC (1), + AC $ )dx 

-b p^ 1 o 


+ AC (1)^ dx - F (y,t) 

-b 1 s 


The Suction force likewise can be separated. The leading-edge inverse- 
square-root pressure singularity strengths for steady and unsteady flow 
can he defined. 


(7-10) C (y) = ^ lim [-/x+b AC (x,y)j 
^ x>-h Po 

(7-11) G.„ (y;t) H lim [- v^x+b AC (x,y;t)] 

^ x^b Pi 

From equations (5-9) and (5-10) , the leading edge suction force in 
terms of the vorticity singularity is 

(7-12) F (y;t) = t p[ 11m (v’xfb Y(x,y;t))]^ 

® a 

where the effect of compressibility is now Included with the /1 -m^ 

a 

factor (deduced from Eq. 12-1 of Ref. 26). Vorticity and pressure dis- 
continuity distributions in the vicinity of the leading edge are related 
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by (Ref. 23, Eq. (5-93)) 


(7-13) ACp(x,y;t) = 


2y 


(x,y;t) 


The suction force in terms of the leading edge pressure singularity 
is, therefore. 


(7-14) F (y;t) = rTp{I [- AC (x,y;t)]}' 

® ^ x-b ^ 


= >^1-Mf np{C (y) + C (y;t))' 
* ^o ^1 


= *^='" np[Cp + 2c c + 4 1 , 

« *^0 *^1 h 


where (7-7), (7-10), and (7-11) have been used. Equations (7-9) and 
(7-14) suggest that chordwlse forces do affect both the steady displace- 
ment solution and the linearized unsteady stability problem. 

Actual computation of the suction force contributions can be 


accomplished through program output of the series coefficients a 


(. 1 ) 

vp 


Insertion of (7-1) into (7-10) and (7-11), for the jth mode, leads to 
the steady and unsteady leading-edge singularity strengths 


(7-15) cj,^\y) = - 2V[’^^ Z f ^^\y) ] 

o V=1 


(7-16) ci^hy) = ’ f^'^^y)]*^ 


v=l 


Here, in taking the limit x -► -b , only the U = 1 chordwise pressure 
distribution terms from (7-4) remain since 
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.(VI) < 


(7-17) lim [/»fb g (x) ] = 
x->^b 


r V* “ 1 

I 0 W > 1 


The single summations that define the modal spanwlse leading-edge 

pressure singularity strength can be arranged for eoniputatlonal purposes 

rr N 
E 

y=l 


(7-18) rcy) = *^3 


(7-19) 


f.(y) = I 


J 


^ V=1 ^ll 


(y) 


This notation, together with summation over all modes, allows the total 
singularity strengths (7-10) and (7-11) to be computed by 


(7-20) C_ <y) = - 2Vv^ I F"(y)q" 
o j ^ ^ 


(7-21) C_ (y,t) = - 2V,^ Z F.(y) q,(t) 

j ^ 


Insertion of (7-20) and (7-21) into (7-14) gives for the suction force 


(7-22) F (y,t) = 4npv"b Tl^ {[!: F’(y)q“r 
s a ^ V V 


+ 2 Z I F"(y)q* q.e^‘^‘^F,(y) 

Vi ^ ^ J 


+ [F F (y) q 
j ^ ■' 


The notation of (7-18) and (7-19) is next adjusted to identify the 
specific structural modes involveid. For the steady problem a rigid 
pitching mode must be used to solve for pressures and loads due to 
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airplane angle of attack, and its associated singularity strength 


parameter is denoted as , For the jth elastic torsion mode F' 


j 


is used; of course, the vertical bending modes Introduce no steady 
lifting loads. Unsteady, leading-edge-slngularity strength parameters 
for the jth elastic torsion and bending inodes are specified by F 
and F . respectively. The computational form of (7-22) is then 


(7-23) F (y;t) = {[F*a + S F® q® Y 

® ^ P=1 % K 


n n 4 -f 

+ 21 [F'ct + Z F“ q? ]F -r*- 
. , a , d) ‘A ^ w, b 

j=l V=1 j 


n n 

. T ^ o 


i«i^t 


+ 2 E [F®a + I ]p. q. e 

j=l “ v=l \ \ 


n w . n * ^ o 

+ [(Z F + I F q )c'^“'']-} 

J=1 ''j ^ .1=1 


The x-component of the ACp--force, as Indicated in (7-9), must also 
be expressed in a form which permits computation. Insertion of the 
modal quantities 

(7-2A) = S f q; 

v=l V 


(7-25) .ACp = ACp ot + E ACp. q“ 

“ V=1 % 


(7-26) 4. = Z f . q. 

‘ .1-1 *1 *1 

Q 

n V . n 

(7-27) ACp, - Z ACp + E ACp 

' W. D 


" 1=1 j 


j = l ^.1 ^ 
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into (7-9) leads to 


(7-28) D (y;t) = 3spV"{[ Z -^v ^Cp“ f, dx « q" 

^ V"1 ^\) 

+ \ ‘t ■>; ‘ + 'A '♦ ■** ^ 

p=l V”! ^p ^\) 1=1 v=l 1 

i=l v-1 ^1 V 

4- E /'l ACp f , dx a q + E E ACp , f dx q“ q. 


n ^ V w* 

+ [ S S / ACp f . dx — ^ q . 
i=ij=i-'’ S *j 

+ E ^ ^a *** 9a 9a 

i=l j -1 *^1 ^1 *1 *.1 - 


where steady, linear unsteady, and nonlinear unsteady terms in the 
generalized displacements are grouped together. 

The chordwlae integrals in (7-28) involving the steady and unsteady 
modal pressures can be computed at any spanwlse station by taking 
advantage of the Rowe program’s capability to compute sectional general- 
ized forces. Program output is of the form 

(7-29) <)Jj(y) = ACp.(x,y)H^(x,y)dx 

where the deflection mode shape H^(x,y) for all integrals in (7-28) 

will be f . (y) . 

*i 

The nonlinear modal equations (4-9) for steady displacements must 
be modified by substitution of the steady compressible lifting airloads 
in place of the strip theory loads (4-3). Inclusion of the steady 
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parts of (7-28) and (7-23) will account for the induced drag caused by 
the three-dimensionality of the flow. Generalized forces which will 
appear in equations (4-9a) and (4-9c) are to be computed directly, but 
the generalized forces for the chordwise bending eauacions (4-9b) 
involve the chordwise forces described above and must be calculated 
separately. These preliminary calculations require the pressure 
coefficients a^^ and sectional generalized forces obtained 

from Rove program output. Program output of generalized forces is of 
the form 


(7-30) ACpj(x.y)H^(x,y) dx dy 


where H^(x,y) is the ith modal displacement and ACp^(x,y) the 
pressure distribution per unit q. . 

Rederlvation of (4-9) with three-dimensional subsonic steadv air- 
loads produces the following nonlinear system: 


(7-31a) 


MPi ‘’w 
a 


MPi 


n n 


TT"?!** 

U-* b 


(t-1) 


[I eh, 

P= 


^ ^"vjU "d 

p = l V=1 ■* ^ 


M=1 v=l i“l p A' v=l ■' V 
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II V 

xMPi \ HPi n n Si 

(7-3it) n‘»; -5^ -ji - h-i) I -f 


n n n v. 

-f- E E ^ Rj 1 — ] 

v=l u=l i-i 


n 


,8. , ,b . „ o 


-2M ^/o f (y)f (y)dy q“ a 

V»1 fi 


(1 < J < n) 


P"1 V=1 '^P V j ^p V 


+ {fhr”)^i dv + 2 z / f;; f® f dv a q? 


V ^ ‘’I > ■ " 

y=l v«x ^p ,1 ^p 

Mi MPi n n n q” q“ 

<73ic, ,■ . (7-.) ^ \ ,. 


ujv,: - £ t„ SS, 


q® q“ n n 


n n T» 


2Ha^'j+2n,o “ -*+2n,\>4-2n ^ 


In (7-31) the notation 2? refars to the steady generalized forces 

j 

computed directly by the Rowe program. Here 2n+l structural modes 
must be used to compute generalized forces: n in vertical bending, n 

in elastic torsion, and one rigid pitching mode. 

Computation of the chordwise load terms in (7-31b) quickly becomes 
an unwieldy task as n is Increased, because of the profusion of 
numerical integrations that become necessary (for n • 3 there would 
be 75 integrations). In order to avoid extensive labor but still retain 
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the most significant effects of the potential chordwise forces, the 

following simplifying assnmptions are adopted. First, integrals 

appearing In the first (j=l) chordwise equation of (7-31b) only are 

retained. Second, only the pressure distribution and slope of the first 

of the elastic torsion inodes are kept; terms containing or f. 

% 

for M > 1 are dropped. Similarly, all terms containing F? for 
U > 1 are removed. This approach should preserve the first-order 
effect of steady chcrdvlse potential forces, yet only five numerical 
spanwlse Integrations will be required. 

Actual computation of the remaining integrals 

^ ^Cp»^(x,y)dxlf^^(y)f^^(y)dy 

is done by direct calculation of the Integrals over % as sectional 
generalized forces with Che Kowe program at eleven spanwlse stations. 

The spanwlse integration is then carried out numerically. For the 
three Integrals 

1 1 
1 1 


the suscaatlons Implied by (7-18]l are first made. This step is followed 
by spanwlse numerical integration. 
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The nc"Linesr solution scheme described In Appendix C was straight- 
forwardly adapted to compute steady deflections in subsonic three- 
dimensional flow from (7-31). Since spanwise induction introduces 
coupling among the torsional modal equations, solution for linear dis- 
placements as an initial estimate becomes an nxn linear matrix problem. 

Equation (4-12) remains as the 3n x 3n linear unsteady modal 
system for stability about the steady equilibrium position. It now 
requires the generalized force matrix to be expressed for three-dlmen- 
sior..-il unsteady compressible flow. As for the steady case, the 
generalized forces relating pressures and displacements in the vertical 
bending and torsion modes can be computed directly by the Rowe program, 
with reduced frequency and Mach number specified. Direct insertion Into 
(4-12a) and (4-12c) is accomplished by the simple substitution, for 
generalized forces and (7^ ^ relating the same two modes. 


(7-32) 


^J.l “ 2TOF -j.i 


Even more so than for the steady airloads, complete inclusion of 
all linear unsteady potential chcrdwlse terms introduces a profusion 
of integrals. Formally, the linear unsteady terms which appear in 
(7“28) and (7-23) enter into generalized forces, computation of which 
Involves spanwise integrations of the terms' products with 
Practically, simplifying assumptions of the type made for the steady 
chordwise terms have been made to keep the number of numerical integra- 
tions at a manageable level. Accordingly, only the first chordwise 
equation of (4-12b) will retain the chordwise force terms. Furthermore, 
only th.e chordwise force integrals containing pressure contributions 
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and displacements for the qj, . q > and q degrees of freedom will 
be kept. (Justification of this approximation is based on the flutter 
behavior observed in Chaptej: VI, which revealed little participation by 
the remaining vertical bending and torsion inodes.) The remaining linear 
unsteady terms from (7-23) and (7-28) are 

(7-33) Dp(y(0 = dx 

H ^ !? h 

+ / ACp f. dx q“ + / , ACp , f, dx q° q (t) 

-h (j)j^ b -b (|)^ 


+ ACp“ f , dx a q . ft) + ACp" f , dx q“ q , (t)} 
-b a (J)^ -b ({)j^ (b^ 


V 


+ SUpV®*^:^ {F® F a -vi + F® F a 

a a b ot W2 b 

'’w 

+ « q. (t) + F? F q? 

q (t) 

W 

+ K Kr < ^ 'I* ^t)} 

<?! «2 *^i ^ *^l ‘^1 *1 


The generalized force matrix terms for inclusion in (4-12b) are then 
found to be 


(7-34) 




/il^ F;(y)F„^(y)f^^(y)dya FJ^ 
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(7-34) 
con’ t 


^l+n,2 


2ii^ ACp^^(K.y)dx]f^^(y)f^^(y)dy 


- F" f dy cc + / F;? F f dy q! } 

k Jl a o a ° *^1 ”’2 '^1 ‘*’1 

‘’l+r.,l+2n = “ 

“D 1 1 


r rb 


8-rrb 


+ r"[/“ACp* dx]E, f dy q? +/„[/“ ACp. dx]f. f dy q” ) 

0 -b ^1 ° n 'I'l "^1 'I’l 


- {/*'■ F" F. f dy a + r F. f dy q” } 

k «• s o a (j)^ o 4 .^ ^c|)^ 


All remaining Q.j,^ . terms in (4-12b) will be xero. 

j+n, 1 

For a given reduced frequency, computation of these generalized 
forces involves program output for both oscillating and steady flow 
conditions. Sectional generalized forces and pressure series coef- 
ficients output by the Rowe program are used, in the same manner 
described above for steady chordwise loads, to compute (7-34). Integrals 
for oscillating flow, of course, are complex. Nine complex spanwlse 
numerical integrations are needed for each k , whereas 126 would have 
been required without the simplifying approximations. Two real integra- 
tions in l+2n iuvolving steady pressures also appear in the 

steady displacement solution. 


B. Flutter Calculation Procedure and Results 

Inclusion of subsonic three-dimensional (3-D) airloads eliminates 
a serious flaw of the strip theory loads used to obtain all Chapter VI 
results. This was the approximate spanwisc load distribution, which is 
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most inaccurate near the tip. Leading-edge suction has been included 
in the same manner as the Chapter V, Section D, analysis of 2-0 flow. 

An effect which was neglected in Che strip theory case, that of the 
x-component of the resultant pressure force normal to the deformed 
chord (cf. Eqs. (5-18) to (5-20) is now retained. It is accounted for 
by the chordwtse terms which are computed with sectional generalized 
forces. The influence of Induced drag upon both steady deformations 
and flutter stability should now be implicitly included by the modeling 
of chordwise loads. Two parameters which must be specified in addition 
to those mentioned in Chapter IV, in order to define a specific wing 

I 

are the aspect ratio (— ) and Mach number. 

b 

The Chapter IV stability calculation method has been modified to 
accept the subsonic 3-D steady and unsteady airloads derived in Section 
A. Since these airloads are now externally generated, iteration of the 
reduced frequency to find neutrally damped eigenvalues is no longer 
feasible. Generalized forces have to be computed beforehand, both for 
steady flow and for oscillatory floxj at preselected reduced frequencies. 
Computation of steady displacements for a given angle of attack and 
flight speed is accomplished as before, now based on Equations (7-31). 
Then for each of the preselected reduced frequencies, the previously 
calculated generalized forces (7-32) and (7-34) are input, the eigen- 
value determinant (4-19) Is assembled, and complex eigenvalues are 
determined. For each of these there are an associated speed and damping 
as shown by (2-25) and (2-27). The speed and damping of each flutter 

mode are then plotted for all k , and a neutrally-damped speed U_ 

r 

is determined by graphical interpolation. The procedure is repeated 
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with a newly estiraated U ; with care, U and can be matched 

e e c 

with Sufficient accuracy after two such steps. 

The expense in computer time to execute the Rowe program and the 
additional effort required to prepare the chordwise loads make It 
desireable to use as few as possible k values. Fortunately, all U 
vs. g interpolation graphs proved to be quite smooth. Flutter speeds 
accurate to three significant digits (sufficient) were reliably obtained. 

The moderate-ar.pect-ratio wing whose stability boundary for 2-D 
airloads appears in Fig. 6-7 has been re-analyzed using 3-D incompres- 
sible aerodynamic theory. Figure 7-1 shows the results,, compared with 
the 2-D flutter calculations. The curve marked "100^ suction" i<7as 
computed with the complete system (4-12), (7-32), and (7-34) whereas 
that marked "Q% suction" was found by repeating the analysis after 
removal of all terms containing the singularity strength parameters 
F, etc. in (7-34). This latter result thus represents the effects of 
forces normal to the deflected airfoil chord only in the dynamic equa- 
tions. The coupling of fore-and-aft bending motions into flutter by 
the leading edge suction forces is therefore absent in the "0% suction" 
case. 


The 2-D stability boundary involves the aerodynamic loads of 
(4-18), which actually do not account for suction. The effect of model- 
ing the suction force, as derived in Section C of Chapter V, is shown 
by Figs. 6-35 and 6-36. These plots suggest that, if a stability 
boundary had been determined with suction accounted for, the 2-0 curve 
in Fig. 7-1 would show less Influence of steady deformation and would 
not drop below U = 7 . The solid, 3-D curve does, of course, account 
for suction. 
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Comparison between the 2-D and 3-D flutter speeds is also made 
more difficult because they involve different steady aerodynamics. In 
order to minimize differences of this kind, the steady tip deflection 
was judged to be the best common measure of steady airload effects. 

Inspection of Fig. 7-1 affirms that the influence of steady 
deformations upon flutter is not Just governed by elastic bending-torsicn 
coupling tut is also sensitive to the manner in which the potential 
aerodynamic loads are applied to the fore-and-aft degree of freedom. 

Even though the chordwise force components represent tilting of a 
relatively largd, approximately vertical resultant force vector, they 
are seen to have a significant stabilizing influence. This conclusion 
follows from comparing the "100% suction" and "0% suction" curves. Any 
analysis of this type will obviously be sensitive to the way in which 
chordwise forces are accounted for. 

The influence of compressibility is next explored by repeating the 
foregoing calculations with airloads computed for Mach numbers of .6 
and .8 . Results are shown in Fig. 7-2. Strictly speaking, this 
procedure involves an Inconsistency, since Mach number Is held fixed 
while velocity is freely varied. The type of calculations required to 
model properly a wing at high subsonic speeds would require iterative 
matching of Mach number and speed for a given flight altitude. 

This refinement is deemed to be excessively costly in both computer 
time and effort. Nevertheless, the reauiLs of Fig. 7-2 provide inter- 
esting qualitative information on Mach number effects. 

As is Increased, the decrease in flutter speed with steady 

deformation for the "0% suction" case becomes less pronounced. The 


221 



results for "100% suction" show less sensitivity to Mach number. In 
either case, no adverse effects due to compressibility upon flutter of 
a lifting wing are revealed, other than those associated with a modest 
decrease In 11^, with increasing . 

When extending this analysis to higher subsonic Mach numbers, it 
must be remembered that the only aerodynamic loads being included are 
those which arise from inviscid, first-order, small perturbation steady 
and unsteady theory for planar lifting surfaces. Induced drag is present, 
as discussed above, but alD. chordwlse forces which arise from either 
direct viscous shears or from modifications to the pressure distribution 
due to the presence of a boundary layer are not modeled. Yet viscous 
effects of this type become increasingly more pronounced as flight 
speeds approach the transonic regime and/or as mean angle of attack is 
increased. 

The chordwise laods associated with direct viscous shear should not 
contain a significant unsteady component, at least .at the lower reduced 
frequencies encountered here. It is expected that they will produce a 
steady drag force aligned with the airfoil section. By contrast, the 
chordwise potential loads represent the horizontal component of a 
relatively large resultant circulatory lift vector. Tilting of this 
vector can introduce unsteady chordwise loads of considerable magnitude. 
Indeed their importance In aeroelastic stability is suggested by 
comparison of the 100% and 0% suction curves in the figures. 

Since viscous shear should cause a predominantly steady drag loading, 
qualitative information regarding its effect on flutter can be inferred 
from the strip theory studies of Chapter VT involving the drag parameter 
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C. For example, no substantial alteration of the dependence of flutter 
speed upon steady lift should be expected from Including this additional 
drag term. The effect of unsteady viscous contributions to aeroelastic 
stability has not been considered anywhere in this investigation. The 
mere prediction of unsteady viscous chordwise loads is still regarded 
as an open question for experimental and analytical research, especially 
when turbulent boundary layers are involved. 

As a final case, the large-aspect-ratio example of Fig. 6-13 has 
been reanalyzed. Aspect ratio is fixed at 20 and Mach number at 0. 
Results are shown in Fig. 7-3. Interestingly, the zero-lift flutter 
speeds for 2-D and 3-D flows are nearly the same; of course, strip 
theory is expected to be more accurate for the larger aspect ratio. 

The difference between the 2-D and 3-D (with suction) stability 
boundaries is actually deceiving, since (as already mentioned) the 2-D 
results do not contain the improved modeling of Che suction force from 
Chapter V. Figures (6-37) and (6-38) suggest the stabilizing effect 
which introducing suction would have upon the 2-D curve in Fig. 7-3. 

The key observation here is that the same type of instability, involving 
substantial participation of the fore-and-aft bending degree of freedom 
observed for 2-D aerodynamic loads, is still observed after 3-D aero- 
dynamics are introduced. Also the Influence of steady deformation upon 
flutter speed remains appreciable. Removal of the unsteady suction 
terms decreases stability, much as was observed for ^ = 10 in Fig. 7-1. 
Points on the stability boundary for oi < .005 could not be found 
reliably because of the extremely light damping in this region; 
approximate curves are shown by dashed lines. 
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In conclusion, the results found in this chapter seem to confirm 
that the phenomena analyzed extensively in Chapter VI are still observed 
when the strip theory airloads are replaced by airloads from 3-D 
aerodynamic theory. Furthermore, the role of chordwise forces due to 
leading edge suction has been found to increase the aeroelastlc stability 
of a wing undergoing steady deformation due to lift. 
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FIGURE 7-3 Flutter Speeds as Influenced by Steady Deformation; Large-Aspect-Ratio Example 
of Fig. 6-13 with = 20 and M = 0 





Chapter VIII 


CONCLUSIONS AND RECOMMENDATIONS 

The following conclusions are drawn from the variety of results 

obtained in this investigation. 

1) The influence of steady drag on flutter speed changes from favorable 
to unfavorable as aspect ratio is increased. The frequency of the 
second transverse bending mode decreases, tending toward the funda- 
mental bending- tors ion flutter frequency, as aspect ratio is 
Increased. 

2) The prediction of the influence of steady drag upon flutter is not 
substantially altered when steady defcrmations due to lift are 
considered. The major effect of steady drag is to reduce divergence 
speed, especially for large aspect ratios. 

3) When a wing has such a large aspect ratio that its fundamental 
fore-and-aft bending frequency is less than the frequency of bending 
torsion flutter at zero steady lift, an instability associated with 
chordwise bending occurs. The critical speeds are lower than the 
zero-lift flutter speed when any steady lifting deformation Is 
present. This type of flutter can disappear at small steady deflec- 
tions when realistic structural damping is introduced, but for 
reasonably large displacements it can still occur. 

A) Steady deformations decrease flutter speed and flutter frequency. 

The effect is most pronounced when the fundamental fore-and-aft 
bending frequency is near the zero-lift bending-torsion flutter 
frequency. 
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5) Thi> aeroclastli: phenomena predicted using Incotopressible strip- 
theory airloads are also observed when three-dimensional, compres- 
sible subsonic airloads are employed. 

6) The inclusion of unsteady leading-edge suction forces moderates 
the predicted decrease in f7iutter speed due to steady lifting 
deformations. 

7) In subsonic compressible flow, with unsteady potential chordvise 
forces, there seem to be only slight adverse effects of H on 
flutter when steady deformations are present. 

One of the predictions of this study is that a high performance 
sailplane undergoing a limit l.oad factor pullout from a dive could 
encounter flutter or divergence at lower speeds than might be antici- 
pated from a conventional aeroelasclc analysis. The increased steady 
drag which accompanies higher would reduce the divergence speed 

considerably, while the deformation of the flexible high-aspect-racio 
wings would change the dynamic aeroelastic stability as well. 

The present analysis could be refined still further. Nonlinear 
aerodynamic effects deserve further attention. Higher mean angles of 
attack would lead to increased importance of the turbulent boundary 
layer, culminating in separation (stall), which will alter the stabiliz- 
ing contribution found due to leading-edge suction in attached flea. 

Vehicles intended to operate at the higher subsonic or low tran- 
sonic speeds (certain RPV's or missiles could still be designed with 
the straight wings considered here) are expected to encounter various 
phenomena which could greatly affect the chorduise loads. Viscous 
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shear, boundary layer modification of the flow, and thickness effects 
with the appearance of shocki^ will all modify the aeroelastic behavior. 

Finally, the transonic and supersonic flow regimes, where drag 
loads are considerably larger than at subcritical speeds, remain largely 
vininvestigated. In these ranges different structural and aerodynamic 
configurations are likely to be involved. 
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Appendix A 


THE ASSUMED MODES 


1. Bending 

The natural mode shapes and frequencies In bending of a uniform 
beam of length £ are found by seeking homogeneous simple harmonic 
solutions of 

El w’"' + nw = 0 

X 

Letting 


the general solution is 


W = A^ sin 


/u 

- y 


+ A„ cos 


— y + A_ sinh 
a 3 


' — y + A, 
a 4 


cosh 



where 



Application of the bending boundary conditions 


vr(0) = W'(0) = W"(£) = M"'(£) = 0 


results in a trancendental equation for the natural frequency eigenvalues 


COS 


cosh — £ 
a 


which are for vertical bending 


(A-i) o> . ^ 



1 = 1.2,3,. 
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are 


The trancendental numbers for 1 _1 i ^ 5 


(A-2) 


= 0.596864162695 
N 2 = 1.494175614274 
- 2,500246946168 
N, = 3.499989319849 
Nj = 4.500000461516 


and the corresponding eignevectors yield the vertical bending natural 
inodes 


(A-3) 



' sin N,n - sinh N.IIl 

— " ■ — 7 rnj — (sinh HN.y - sin IlN.y) 

cosh TIN. + cos HN^ I I-' I-" 

1 i ' 


+ cosh riN.y - cos nN.yl 
1 i 


expressed In orthonormal form so that 


1 


1 


The modal property in the steady equations (4-9) is related to 

the modal integral 


(A-4) 


1 

/ f dy 
o 


.. PslnirN. - sinh srN, 
^ I 1 i 

Nj^IT l^cosh irN^ + cos 



Fore-and-aft bending natural modes f are the same as in (A-3), 

^’i 

however the natural frequencies are 
(A-5) 
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The f are Illustrated in Fig. A-1 for 1 < 1 < 5 . 
w — — 

1 

2. Torsion 

Natural mode shapes and frequencies in torsion are simply 
(A-6) £. » sin II(j-Si)y 

1 

(A-7) 

vfhich result from the elementary Sturm-Liouvllle problem for torsion 
of a uniform rod. These modes are not normalized since 

f? (v)dv = h 
0 - 
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Appendix 8 

CALCULATION OF MODIFIED BESSEL FUNCTIONS 

The modified Bessel functions Kj^(s) appearing in 

the expression (5-1) for the generalized Theodorsen function can be 
computed by using the following ascending power series expansions, 
drawn from Abramowitz and Stegun (Ref. 26, Equations 9.6.10, 9.6.11, 
and 9.6.13). 


(B-1) 


Iy(i) = (3is)^ 


7 

j!r(v+J+l) 


(B-2) 


K (s) = - {En(Jss) + Y.}l (s) 

O 6 0 


I I ni'o I ni' I 

OTT ~2ly^ (H-5+r)-(3jjT- + 


(B-3) K^(i) = !s(Jii)“^ + £ndi)I^(i) 


_ ” /ur2^j 

- k(hs) 1. {ij^(j + 1) + ij)(j + 2)} ^ ^ j , 


where 


r(v + 1) = v! 


l(i(l) =-y 

-1 

'f'(v) = -y + IJ , v>2 

j=l 

Y = 0.5772156649 (euler's constant) 

e 


No convergence difficulties with the power series expansions were 


encountered over the range of s that occurred In this investigation. 
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Simple harmonic 
frequency C(k) 
substitution s 


alrloac’b involving the Theordorsen function of reduced 
were computed by the same procedure, simply by the 
= 0 + ik . 
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Appendix C 

SOLUTION FOR NONLINEAR STEADY DISPLACEMENTS 


With the vector of generalized disnlacements defined as 


(C-1) 


q. 


tq: 


n 1 n 


’q* ^ 

n 


the nonlinear equations (4-9) can be expressed. 


(C-2) F(q“) = 0 


where 

(C-3) 


F = {F 


.F ,F , 

w ’ V ’ 

n 1 


V 0 
n 


, ...,F^ } 


Let an initial estimate q? <l be found by solving the linear steady 

equations (4-20), then linearize F about q, » by first-order Tavlor 

■w(o) 

series, 

(C-4) - q;„)> 4- H.O.T. 

The Jacobian matrix . J contains partial derivatives of the F's with 
respect to the q°'s evaluated at i(o) ; its elements are shown in 
Fig. C-1. Equation (C-4) can be used to solve for (q/v, - q/n\5 by 
the linear approximation 

(c-5) * '"Wo)>’<aa) -i(c)> -o 

and the first iteration solution is 


(C-6) 





) + q 

/V"' 


(o) 
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This process is repeated until satisfactory convergence Is achieved, 
with the general jth step given by 


(C-7) 




( 3 )^ 




L(j+i) 


(q 


0 + 1 ) 




O) 


In practice, more Chan four or five Iterations rarely were required, 
for cases near static divergence speeds. 


2A0 



w MPi MPl n n 

^1,1 ■ bii^ «ij <"-» Eo^ ^ ^ "ijw i; 

J' -> J y=X V»1 ^\) 


3F 


w, MPl 


j,i+Ji aq* 




3F 


MPl n 


n n 


^J,l+2n - 3,- * - <’-W »11„ b“ - “ “vjUl b'’’ - 2'lJ 


a? 


j+n,i aq^ 


KPi 

0 


(T-1) Z 


vfi \ 


3F 


MPl. 


j+n, 1+ u aq. 


A- = xn^N*: 6 ' ® 


MPl. n n 


j bU^ 11 


(T-1) 


z R. . q! q! 
bU , , liuv d> 

v=l y=l 


ar 


"V MPl n ‘iw n n 'v 

^j+n,i+2n ' aq»^ ^ ‘ ■*■ ^ \ 

^^<() n n '^v n 

^j+2ti,l+n = 3q» “ f" '**V ~ 


3F 


<t>^ „ Mi 

i = ~ A) S 


■^j+2n,i+2n " 9q® U 


ij 


MPl^ n n 


W " " 


q o ^ O 

q 

V 

u V, 


+ (T-1) [ E E a E E R,^^i-h^-j- 

" v=l y=l ^ ^ v=l u=l ^ 


ij 


,1 (l=j) 

lO (l3<j) 


(1 < 1 < n) , (1 < i < n) 


FIGURE C-1 Elements of the Jacobian Matrix 
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